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We present a detailed theoretical overview of the thermodynamic properties of the dipolar spin 
ice model, which has been shown to be an excellent quantitative descriptor of the Ising pyrochlore 
materials Dy2Ti207 and Ho2Ti207. We show that the dipolar spin ice model can reproduce an 
effective quasi macroscopically degenerate ground state and spin-ice behavior of these materials 
when the long-range nature of dipole-dipole interaction is handled carefully using Ewald summation 
techniques. This degeneracy is, however, ultimately lifted at low temperature. The long-range 
ordered state is identified via mean field theory and Monte Carlo simulation techniques. Finally, we 
investigate the behavior of the dipolar spin ice model in an applied magnetic field, and compare our 
predictions with experimental results. We find that a number of different long-range ordered states 
are favored by the model depending on field direction. 



I. INTRODUCTION 

A. Water Ice and Spin Ice 

Frustrated or competing interactions are a common 
feature of many condensed matter systems^ In mag- 
netic materials, frustration arises when the system can- 
not minimize its total classical ground state energy 
by minimizing the energy of each spin-spin interaction 
individuallyi2i2i4 When competing interactions cannot be 
simultaneously satisfied as a consequence of the arrange- 
ment of spins on a geometrical unit, such as a triangle 
or a tetrahedron, a system made of an assembly of such 
units is said to be geometrically frustrated. Geometric 
frustration has been studied extensively in recent years, 
with the discovery of classical systems that do not dis- 
play any ordering or dynamical phase transitions down to 
the lowest temperatures (for recent reviews see Refs. |5|- 
ITll) . Furthermore, much current research effort is being 
deployed to investigate the exotic behavior of quantum 
frustrated systemsi 12 i 13 i 14 i 15 In highly frustrated systems, 
weak quantum fluctuations may work to select a unique 
ground state that is not stabilized at the classical level, 
while strong quantum fluctuations (e.g. small spin num- 
ber value, S) can give rise to novel quantum disordered 
states^ Real materialiLiLiS and model systems with 
strongly correlated electrons in the presence of strong 
magnetic frustration display interesting exotic properties. 

While geometric frustration most commonly arises be- 
tween spins interacting antiferromagnetically (AF), Har- 
ris and collaborators^^ 1 " showed that the pyrochlore lat- 
tice of corner sharing tetrahedra with Ising spins point- 
ing along a local cubic (111) axis constitutes a new class 
of geometrical frustration when nearest neighbor inter- 
actions are ferromagnetic (FM) (See Fig. ^jS 2 . 3 As a 
consequence of the frustration on this lattice, the Ising 
pyrochlore ferromagnet has a lowest energy ground state 
configuration that is very closely analogous to an entirely 
different yet very common frustrated condensed matter 




FIG. 1: The (fff) Ising pyrochlore lattice. The lower left 
"downward" tetrahedron of the pyrochlore lattice shows Ising 
spins as arrows. Each spin axis is along the local (111) quan- 
tization axis, which goes from one site to the middle of the 
opposing triangular face (as shown by the disks) and meets 
with the three other (111) axes in the middle of the tetrahe- 
dron. For clarity, black and white circles on the lattice points 
denote other spins. White represents a spin pointing into a 
downward tetrahedron while black is the opposite. The entire 
lattice is shown in an ice-rules state (two black and two white 
sites for every tetrahedron). The hexagon (thick gray line) 
shows a minimal size loop move, which corresponds to revers- 
ing all colors (spins) on the loop to produce a new ice-rules 
state. 



system — namely water iceAii In the low temperature 
— low pressure phase of water ice (the so-called "hexag- 
onal ice" , phase Ih), the oxygen atoms are arranged on a 
hexagonal lattice, each oxygen having four nearest neigh- 
bors. Bernal and Fowler 24 and Pauling 25 were the first 
to propose that the hydrogen atoms (protons) within the 
H2O lattice are not arranged periodically, but are disor- 
dered. These hydrogen atoms on the O—O bonds are not 
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FIG. 2: The local proton arrangement in ice, showing oxygen 
atoms (large white circles) and hydrogen atoms (small black 
circles) arranged to obey the ice rules. The displacement of 
the hydrogen atoms from the mid-points of the oxygen-oxygen 
bonds are represented as arrows, which translate into spins on 
the pyrochlore lattice in Fig. Q 

positioned at the mid point between both oxygen atoms, 
but rather each proton is (covalently) bonded "near" one 
oxygen and (hydrogen-bonded) "far" from the other such 
that the water solid consists of hydrogen-bonded H2O 
molecules (see Fig. EJ. In the Pauling model, ice Ih is 
established when the whole system is arranged according 
to the two ice rules: 

1. Precisely one hydrogen atom is on each proton 
bond that links two nearest neighbor oxygen atoms. 

2. Precisely two hydrogen atoms are near each oxy- 
gen atom (spin in) and two are far (spin out) (see 
Fig. |3 . 

A consequence of this structure, and the subsequent ice 
rules, is that there is no single unique lowest energy state. 
Indeed, there exists an infinitely large number of degen- 
erate low energy states that fulfill the ice rules and, if 
the degeneracy was truly exact, would manifest itself 
as a residual entropy at zero temperature (called zero 
point entropy). Linus Pauling 2 ^ estimated theoretically 
the residual entropy, S(T — > 0), of ice as 

<?(T^0)«|ln^ , (1) 

where R w 8.31 J mole _1 K _1 is the molar gas constant. 
Pauling's result is not exact, but is accurate to within a 
few percent compared to experiments^ 

Returning to the magnetic Ising pyrochlores, the anal- 
ogy to water ice arises if the spins are chosen to repre- 
sent hydrogen displacements from the mid-points of the 
0—0 bonds (Fig.[2J). The ice rules of two protons close, 
two protons further away corresponds to the two spins 
in — two spins out configuration of each tetrahedron on 
the pyrochlore lattice. Because of this direct analogy 
between water ice and the Ising pyrochlores, Harris et 
ali 20 i 21 called the latter spin zcei 9 i n i 22 i 23 We note, how- 
erver, that common water ice at atmospheric pressure, 
ice Ih, has a hexagonal structure while here, the magnetic 
lattice has cubic symmetry. Strictly speaking, the Ising 
pyrochlore problem is equivalent to cubic ice, and not 



the hexagonal phase. Yet, this does not modify the "ice- 
rule" analogy (or mapping) or the connection between 
the statistical mechanics of the local proton coordination 
in water ice and the low temperature spin structure of 
the spin ice materials. 

An important point must be emphasized here. In 
both ice water and spin ice, the microscopic origin of 
the residual zero point entropy arises from the "simplic- 
ity" and "under-constraints" in the problem. Indeed, 
the constraints (rules) to construct a minimum energy 
ground state, which arise from the underlying micro- 
scopic Hamiltonian, are so "simple" that an infinite num- 
ber of configurations of the dynamical variables at stake 
(proton position in ice, and spin direction in spin ice) 
can be used to make a minimum energy state from which 
the extensive residual ground state entropy S(T — ► 0) 
results. 



B. Dipolar Spin Ice 

Experimentally, it is known that the single ion ground 
states of the rare earth ions Dy 3+ and Ho 3+ in the 
pyrochlore structure are described by an effective clas- 
sical Ising doublet i - ° ,2 i Specific heat measurements by 
Ramirez 29 on the compound Dy2Ti20y have shown that 
the "missing" magnetic entropy not recovered upon 
warming the system from T sw 0.4 K to 10 K, agrees 
reasonably well with Pauling's entropy calculation above, 
S ~ S(T — > 0), thereby providing compelling thermody- 
namic evidence that Dy2Ti20y is a spin ice material 30 
(see Fig. [3| . While early neutron scattering and mag- 
netization measurements first suggested that Ho2Ti20y 
was a spin ice material^ some subsequent specific heat 
measurements and numerical simulations by Siddharthan 
and co-workers were interpreted as evidence for a freez- 
ing transition to a partially ordered state as opposed 
to spin ice behavior in that materialiii* 3 ^ 3 ^ However, 
more recent specific heat magnetization 3 ^^ and 
neutron scattering experiments^ supported by Monte 
Carlo (MC) simulations, 33 appear to confirm the ini- 
tial proposal^ , that H02T12O7 is indeed a spin ice ma- 
terial akin to Dy2Ti2 0y. Other magnetization measure- 
ments have recently been reported that also argue for 
spin ice behavior in the closely related Ho2Sn2 Orr^^ 
Dy2Sn207, 3S and H02RU2O7 39 materials. The dynami- 
cal properties of these materials at the spin ice freezing 
point appear somewhat puzzling and are the subject of 
an increasing number of studies4£L£i42i22i 

Following the initial spin ice proposal in 1997 by Harris 
and co-workers )22i2i it appeared that the spin ice mate- 
rials obeyed the simple ferromagnetic nearest neighbor 
model mentioned above. This model intuitively gives 
rise to a degenerate spin ice ground state because of the 
equivalent energies of the six different tetrahedron config- 
urations that make up the ground state of this geometri- 
cally frustrated unit. However, the nearest neighbor spin 
ice model is too simple to accurately describe the physical 
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FIG. 3: (a) Specific heat and (b) entropy data for Dy2Ti207 
from Ref. l29t compared with Monte Carlo simulation results 
for the dipolar spin ice model, with J nn = — 1.24K and D nn = 
2.35K. 



properties of real materials composed of the rare-earth 
ions Ho 3+ and Dy 3+ (see Ref. l3lf) . Firstly, the mag- 
netic cations Ho 3+ and Dy 3+ in Ho 2 Ti 2 C"7 and Dy 2 Ti 2 07 
carry a large magnetic momentj22i2£ [i, of approximately 
IO^b- This entails strong magnetic dipole-dipole interac- 
tions in these materials. Indeed, the strength of the dipo- 
lar interaction at nearest neighbor distances, D n n, is of 
order 2 K, which is of the same order of magnitude as the 
overall magnetic interaction energy scale in these materi- 
als as estimated by the Curie- Weis temperature, 9qw ~ 1 
K, extracted from DC magnetization measurements. Sec- 
ondly, rare-earth ions possess very small exchange ener- 
gies, which is roughly the same order of magnitude as 
0cw and D nn . Consequently, dipole-dipole interactions 
in Ho 2 M 2 7 and Dy 2 M 2 7 (M=Ti, Sn) are very sig- 
nificant and constitute an order one energy scale in the 
problem. This is the reverse of what is observed in tran- 
sition metal compounds, where the exchange interaction 
predominates and the dipolar interaction can be treated 
as a very weak perturbation. Finally, the nearest neigh- 
bor exchange interaction in Ho 2 Ti 2 07 and Dy 2 Ti 2 07 is 
actually antifcrromagnetic, which would by itself cause 
a phase transition to a Neel long-range ordered q = 
stat e - 21 ' 23 (see Fig. |SJ). Consequently, we consider the 
simplest model of (111) Ising pyrochlore magnets with 
both nearest-neighbor exchange and long-range magnetic 
dipole-dipole interactions with the Hamiltonian: 
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Here the spin vector S° = afz a labels the Ising moment 
of magnitude |S"| = 1 at FCC lattice site Ri and tetrahe- 
dral sub-lattice site coordinate r a , where the local Ising 
axis is denoted by z a and the Ising variable is af = ±1. 



The vector R: 
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r connects spins S£ and S^. J 



represents the exchange energy and D the dipolar energy 
scale (J > and D = in the spin ice model originally 
proposed by Harris el alm^ which we refer to as the "near 
neighbor spin ice model"). Because of the relative local 
(111) Ising orientations, the nearest neighbor exchange 
energy between two spins is J nn = J/3. The dipole in- 
teraction is calculated from 
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Experimentally, from magnetization measurements^! and 
analysis of the crystal-field levels via inelastic neutron 
scattering, 28 it is known that the moments of the Dy 3+ 
and Ho 3+ rare-earth ions in the pyrochlore lattice are 
fi w 10/iB, and the nearest neighbor distance r nn is 
approximately 3.54 Angstroms. From Eq. we get 
the dipole-dipole interaction at nearest neighbor dis- 
= 5D/3, since z a ■ i b = -1/3 and 
- , = -2/3 in Eq. ©. For both Ho 2 Ti 2 7 
and Dy 2 Ti 2 7 , Am ~ 2.35K. 

In order to consider the combined role of exchange and 
dipole-dipole interactions, it is useful to define an effec- 
tive nearest neighbor energy scale, J c fr, for (111) Ising 
spins: 
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where J nn = J/3 is the nearest neighbor exchange energy 
between (111) Ising moments. This simple near-neighbor 
description of the system suggests that a (111) Ising sys- 
tem could display spin ice properties, even for antiferro- 
magnetic nearest neighbor exchange, J nn < 0, so long as 
Jeff = Jnn + -Dnn > 0. Fits to experimental data give 

J nn 0.52 K for Ho 2 Ti 2 7 & and J nn 1.24 K for 

Dy2Ti 2 07i^ Thus, J e ff is positive (using D nn = 2.35K), 
hence ferromagnetic and frustrated, for both Ho 2 Ti 2 07 
(J G ff ~ 1.8 K) and Dy 2 Ti 2 7 (J cff - 1.1 K). It would 
therefore appear natural to ascribe the spin ice behav- 
ior in both Ho 2 Ti 2 C>7 and Dy 2 Ti 2 C>7 to the positive J e ff 
value as in the simple model of Bramwell and Harris. 21 
However, the situation is more complex than it appears. 

Dipole-dipole interactions are "complicated": (i) they 
are strongly anisotropic since they couple the spin, S°, 
and space, R^ b , directions, and (ii) they are also very 
long-ranged (oc IR^I" 3 ). For example, the second near- 
est neighbor distance is \/3 times larger than the nearest 
neighbor distance, which means that the second near- 



est neighbor dipolar energy is D n 



0.2D„ 



This 
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implies an important perturbation compared to J e s = 
^nn + Am < -Dim j especially for antiferromagnetic (nega- 
tive) J nn . Specifically, for Dy 2 Ti 2 07, the second nearest 
neighbor energy scale is about 40% of the effective near- 
est neighbor energy scale, J c g , a large proportion! There- 
fore, one might have expected that the dipolar interac- 
tions beyond nearest neighbor would cause the different 
ice-rules states to have different energies, hence possibly 
breaking the degeneracy of the spin ice manifold, similar 
to what happens in the kagome^ and pyrochlore Heisen- 
berg antiferromagnetsa& when exchange interactions be- 
yond nearest-neighbor are considered. In Eq.(|2J, if the 
dipolar term is summed beyond nearest neighbor, one 
could expect a long-ranged Neel ordered state at a crit- 
ical temperature Tn ~ 0(D nn ). Thus, here arises one 
of the main puzzling and interesting problems posed by 
the dipolar spin ice materials that can be summarized by 
two questions: 

1. Are the experimental observations of spin ice be- 
havior in real materials consistent with dominant 
long-range dipolar interactions? 

2. If so, why do long-range dipolar interactions fail to 
destroy spin ice behavior and give rise to long-range 
order at a temperature Tn ~ 0(D nn ) ? 

Results from Monte Carlo simulations on the dipolar 
spin ice model attempting to answer the first question 
above were first reported in Ref. and Ref. |U In 
that work, the dipole-dipole interactions were cut-off at a 
distance of five^i or ten and twelve nearest-neighbors. 32 
In those studies the thermodynamic behavior was found 
to be consistent with spin ice behavior for a model of 
Dy2Ti2 07, provided the exchange interaction was made 
to extend far beyond nearest neighbor^ but not for a 
model of Ho 2 Ti207. A subsequent work^4 considered 
the Hamiltonian of Eq. © with only nearest-neighbor 
exchange and the value of J as an adjustable parame- 
ter. In that work, the long-range dipole-dipole interac- 
tion was handled using the well-known Ewald method, 
which derives an effective dipole-dipole interaction be- 
tween spins within the cubic simulation cell. The Monte 
Carlo simulations were carried out by slowly cooling the 
simulated lattice, subject to the usual Metropolis algo- 
rithm. Numerical integration of the specific heat di- 
vided by temperature was performed to determine the 
entropy of the system44 For a parameter J appropri- 
ate for the Dy 2 Ti 2 07 spin ice material (see Section II 
below), the dipolar spin ice model retained Pauling's en- 
tropy (Eq. |JT]|). in good agreement with experiments on 
h^TbCv (Fig. |3J). Following the same approach as in 
Ref. |4J, recent Monte Carlo simulations have found good 
agreement between the dipolar spin ice model, specific 
heat measurements and elastic neutron scattering, as well 
as experiments on Ho2Ti207. 33 Finally, mean-field the- 
ory calculations of the neutron scattering intensity, valid 
in the (paramagnetic) temperature regime T 3> #cw that 
consider large distance cut-off of the dipole-dipole inter- 
actions have been found to be in good agreement with 



experiments on H^S^O? 3 ^ and Ho2Ti207;— Conse- 
quently, there is now strong compelling evidence that the 
long-range dipolar interaction is responsible for the ice 
behavior and the subsequent retention of zero-point en- 
tropy in rare-earth based insulating pyrochlore magnets* 



C. True Long Range Order at Low Temperature in 
the Dipolar Spin Ice Model 

Having answered question Jfl above in the affirma- 
tive, one is then faced with addressing question #2. 
The Monte Carlo results mentioned above^ show that 
spin ice behavior arises from the combination of nearest- 
neighbor exchange, J lm , and dipole energies, D nn , which 
create an effective ferromagnetic Ising model J c g at near- 
est neighbor as long as J e g = J nn + D nn > 0, akin 
to Harris and BramwelPs simple nearest-neighbor Ising 
mO( j e l i 20 4 2i i 49 4 50 However, the long-range dipolar interac- 
tion does not appear to destroy the spin ice degeneracy 
(and subsequent retention of zero point entropy) created 
by this effective ferromagnetic nearest neighbor inter- 
action. In support of this picture, a mean-field theory 
(MFT) calculation finds that remaining (beyond nearest 
neighbor) dipole-dipole interaction terms, which couple 
every spin in the system with varying strength depend- 
ing on their separation distance, is "screened" to a large 
degree^i*^ This means that the degeneracy between dif- 
ferent ice-rules obeying states is almost exactly fulfilled 
by carefully including the long distance dependence of the 
dipolar term in the Hamiltonian. However, and perhaps 
most interestingly for these Ising pyrochlore systems, the 
same mean-field calculation suggests that the screening 
of the long-range terms is not perfect, and that the asso- 
ciated spin ice manifold is only guasi-degenerate, due to 
some small remaining effective energy scale, and that a 
unique ordering wavevector is selected This sug- 

gests that at some temperature below the onset temper- 
ature of spin ice correlations, the dipolar spin ice model 
should in principle favor the unique long-range ordered 
state selected by this remaining ("unscreened") pertur- 
bative dipole-dipole energy. 

One might naively expect that such an ordered state 
should be found in the MC simulations However, 
this does not happen, as measurements of the tempera- 
ture dependent acceptance rate of the simulations make 
it apparent that the standard single (Ising) spin flip 
Metropolis algorithm experiences a dynamical "freezing" 
at a temperature ~ 0.3 K for J nn and D nn parameters ap- 
propriate for Dy 2 Ti 2 07ii and T w 0.6 K for Ho 2 Ti 2 7 i2a 
If the dipolar interactions are cut-off at some arbitrary 
distance, R c , one can generate scenarios where, depend- 
ing on specific numerical values for J nn ,-D nn and R c , a 
selected state is dynamically accessible before the spin- 
ice manifold freeze-out, as was found in simulations where 
dipole interactions are cut-offi 3 ^* 3 ^ Consequently, akin to 
the approaches used in ice lattice models^SiS one must 
introduce non-local dynamics in the simulation to com- 
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FIG. 4: The long-range ordered q = (0, 0, 2ir /a) dipolar spin 
ice ground state. Projected down the z axis (a), the four 
tetrahedra making up the cubic unit cell appear as dark gray 
squares. The light gray square in the middle does not repre- 
sent a tetrahedron, but its diagonally opposing spins are in 
the same lattice plane. The component of each spin parallel 
to the z axis is indicated by a + and - sign. In perspective (b), 
the four tetrahedra of the unit cell are numbered to enable 
comparison with (a). 



bat this freezing-out and maintain simulation equilibrium 
down to lower temperatures. The inclusion of non-local 
"loop moves" in the dipolar spin ice model promotes the 
development of a long-range ordered phase via a sharp 
first order phase transition at T w 0.18 a much 

lower temperature than the onset temperature for spin 
ice correlations at T ~ 1.2 K in Dy 2 Ti 2 Q , 2 ?i 44 and 
T ~ 1.9 K in Ho 2 Ti 2 7 . 33 The ground state found in 
the loop MC simulations has zero total (bulk) magneti- 
zation (recall that each tetrahedron individually carries 
a net magnetic moment in each of the ice-rule obeying 
states). See Fig. 01 for the spin configurations in this 
ground state. The pre-transitional specific heat and the 
latent heat associated with the first order transition re- 
covers all of Pauling's missing entropy in the model. The 
ordered state that is found in the loop MC simulations 58 
corresponds to the ordered state predicted by mean field 
theoryi 5 ^ 5 ^ In other words, the dipolar spin ice model 
possesses on its own, without invoking energetic per- 
turbations and/or thermal and quantum fluctuations, a 
unique (up to trivial global symmetry relations) classical 
ground state with zero extensive entropy. 

Also, using MC simulations and direct Ewald energy 
calculations, we investigate the behavior of the dipolar 



spin ice model in an external magnetic field. With appli- 
cation of a large field along three different crystal symme- 
try directions, three different long-range ordered ground 
states appear. With large fields parallel to the [100] crys- 
tal direction, the ground state is the ice-rules q = 
structure identified by HarrisiSI For large fields parallel to 
[110], the ground state is the ice-rules q = X stated and 
for large fields along [111], the ice rules are broken and 
a three-spin in, one-spin out spin configuration becomes 
the lowest energy state. The experimentally determined 
field dependence of the magnetization and specific heat 
for fields along the [100], [110] and [111] directions in 
the Dy 2 Ti 2 C>7 spin ice material agree quantitatively well 
with the Monte Carlo results for the long-range dipolar 
spin ice modeh fin ' fi1 i 62 ' fi3 i 64 



D. Phases of Dipolar Spin Ice 

Using Monte Carlo simulations, the phase diagram for 
the dipolar spin ice model can be mapped out (Fig. [SJ. 
To summarize the results, spin ice correlations develop 
for all cases where the effective nearest neighbor energy 
scale Jeff /-Dim > 0.095 (ferromagnetic), and the tem- 
perature is below the broad peak in the specific heat, 
Tpcak- For T/D nn < 0.08, independent of the value of 
J nn (as long as J e g/D nn > 0.095), the system orders 
into the long-range ordered state, with the help of the 
loop moves in the simulation. For J nn /D nn less than 
-0.905 (J e g/D Dn < 0.095), the system orders into an an- 
tiferromagnetic q = Neel ground state, where every 
tetrahedron in the system has an all-in or all-out spin 
configuration at low temperaturesiSiiS 3 ^^^ The region 
around J n n/D nn = — 1 shows hysteresis at low tempera- 
tures. Because of the close cancellation of energy scales, 
we imagine that real materials which fall into this region, 
e.g., Tb 2 Ti 2 7) 67 i 68 ' 69 i 70 will be particularly susceptible 
to the influence of small perturbations (such as exchange 
beyond nearest-neighbor or finite, as opposed to infinite, 
Ising anisotropjif^i) with the result of possible ordering 
into long-range ordered states 7 ^ distinct from the two 
shown in Fig. 



E. Outline 

The rest of the paper is organized as follows. In the 
next section we present results from conventional single 
spin flip Monte Carlo simulations that show how spin ice 
behavior develops at finite temperatures in the dipolar 
spin ice system whenever the effective nearest-neighbor 
coupling is ferromagnetic (J e g = J n n + -D nn > 0). Results 
from mean-field theory are presented in Section III. There 
we show that there exists a weak selection of a unique 
ordering (critical or soft mode) at q = (0,0,27r/a). Mo- 
tivated by our mean-field results, we undertake a numer- 
ical search for a long-range ordered state in the model of 
Eq. J2J). Section IV discusses the details of a loop Monte 
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FIG. 5: The phase diagram for the dipolar spin ice model. 
The antiferromagnetic ground state is an all-spins-in or all- 
spins-out configuration for each tetrahedron. The spin ice 
configuration, which includes the q = (0,0, 2n/a) ground 
state, is a two spins in-two spins out configuration for each 
tetrahedron. The region encompassed between the quasi ver- 
tical dotted lines displays hysteresis in the long-range ordered 
state selected (q = vs. q = (0, 0, 2n/a)) as Jnn/Dnn is var- 
ied at fixed temperature T. 



Carlo algorithm that avoids the freezing phenomenon ob- 
served in a MC simulation employing local single spin flip 
dynamics. Section V presents the detailed results from 
the loop Monte Carlo simulations. The results for the 
field dependence of the ground state energy and mag- 
netization for fields along [100], [110] and [111] are pre- 
sented in Section VI. We conclude the paper with a brief 
discussion in Section VII. We have included several ap- 
pendixes. Appendix 1^1 contains a discussion of the Ewald 
technique for dipolar interactions in real space (MC sim- 
ulations) and in momentum space (MFT). In Appendix 
IbI the q— dependent susceptibility, to quadratic order, is 
derived via high temperature series expansion to demon- 
strate the connection of the mean-field T c to the onset of 
long-range correlations. Appendix El discusses some of 
the effects of a finite demagnetization factor on specific 
heat results. 



simulation cell. 73 Unlike in dipolar fluid simulations, 74 
the pyrochlore lattice constrains the positions of the spins 
in the simulation, allowing the Ewald interactions to be 
calculated only once, after which a numerical simulation 
can proceed as normal. Appendix lAl contains a brief dis- 
cussion of the Ewald method applied to MC real space 
simulations. 

Simulations on the dipolar spin ice model are carried 
out using the standard single spin flip Metropolis algo- 
rithm. To mimic the experimental conditions pertinent 
to real materials, the simulation sample is cooled slowly. 
At each temperature step, the system is equilibrated 
carefully, then thermodynamic quantities of interest are 
calculated. Since D nn can be determined once the crystal 
field structure of the magnetic ion is known, the nearest 
neighbor exchange J nn is the only adjustable parameter 
in the model. The single spin flip Metropolis algorithm 
is able to map out three different regions of the phase 
diagram shown in Fig. [S] Thermodynamic data indi- 
cates that when the nearest neighbor exchange is AF and 
sufficiently large compared to the dipolar interactions 
(J nn < and | J nn | ^> -D n n), the system undergoes a sec- 
ond order phase transition (in the three dimensional Ising 
universality class) to an all-in or all-out q = ground 
state. In the spin ice regime, J c g = J nn + D nn > 0, 
each specific heat data set for different J nn shows quali- 
tatively the same broad peak as observed in the nearest 
neighbor FM exchange model^ which vanishes at high 
and low temperatures (see Fig. EJ. The height, C pea k, 
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II. THE DIPOLAR SPIN ICE MODEL: 
CONVENTIONAL METROPOLIS MONTE 
CARLO 

In this Section we present the results of Monte Carlo 
simulations of the dipolar spin ice Hamiltonian [21 using 
a standard single spin flip Metropolis algorithm. To use 
Eq. (J2J within a simulation, the dipole-dipole interac- 
tion must be handled with care. A lattice summation 
of such interactions is conditionally convergent due to 
its 1/R 3 nature. In order to properly handle the long- 
range nature of this term, we implement the well known 
Ewald method in the simulations, which derives an effec- 
tive dipole-dipole interaction between spins within the 



FIG. 6: Specific heat for system size L — 2, with tempera- 
ture, T, re-scaled into units of the effective nearest neighbor 
interaction J e g = J nn + D nn . J n n/D nn = corresponds to 
purely dipolar interactions, while J n n/D nn = oo corresponds 
to nearest neighbor FM exchange only. Simulation runs for 
L — 4 were also performed, but revealed no important finite 
size effects. 

and peak temperature position of this peak, Xp 0a k, show 
very little dependence on system size, for simulation cells 
of L = 2,3,4,5 and 6. However, both C pea k and T pca k 
are found to depend strongly on the ratio of J nn /D nn , as 
illustrated in Fig. [7| 

By re-scaling the temperature scale for the specific heat 
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FIG. 7: Dependence of the simulated specific heat peak height 
Cp C ak and temperature location of C poa k and T pca k on ex- 
change and dipole-dipole interaction parameters. In this fig- 
ure D nn is set to 2.35K. 



corresponding to a number of different interaction param- 
eters J nn and -Dnn j one can expose more clearly the de- 
pendence of the specific heat on the competition between 
the nearest neighbor exchange J nn and the dipole-dipole 
interactions. In Fig. this dependence is illustrated in 
the regime J nn / D nn > 0. This figure shows that in terms 
of an effective energy scale, J c s, the medium to long- 
range effects of the dipolar interactions are in some sense 
"screened" by the system, and one recovers qualitatively 
the short range physics of the nearest neighbor spin ice 
model. As the nearest neighbor exchange interaction be- 
comes AF (see Fig. [7] for J nn /D nn < 0), we find that the 
approximate collapse onto a single energy scale becomes 
less accurate, with the specific heat becoming dependent 
on J nn / D nn . It is within this regime that we believe that 
both Ho2Ti207 and Dy2Ti20y are realized, as we now 
discuss. 

Since D nn is calculated from Eq. iJSJ, J nn must be de- 
termined from experimental data. By fitting either the 
height Cp 0a k or the peak temperature T peak of the maxi- 
mum of the specific heat curves of the Monte Carlo sim- 
ulation to the experimental results^ (Fig. [3^,), we ob- 
tain a value of J nn = — 1.24K for Dy2Ti2C>7. The re- 
sults of this fitting are illustrated in the top panel of 
Fig. [3J A fitting of the height or peak temperature of 
the experimental magnetic contribution to specific heat 
for Bx>2Ti207 gives J nn = — 0.52K for this material^ 
Contrary to what is reported in Ref. |3^, we, therefore, 
conclude that Ho2Ti 2 07 is "deeper" (J e g more positive 
for Ho2Ti207 than for Dy2Ti2C>7) in the spin ice regime 
(farther to the right in Fig. than Dy2Ti2C>7 . As ini- 
tially reported in Ref. |7^, the temperature dependence 
of the specific heat for Bx>2Ti2 07, is less straightforward 
to interpret than for Dy2Ti207i2 9 " In Bx>2Ti207 the spe- 
cific heat possesses an important contribution from a nu- 
clear component due to a large hyperfine splitting of the 
nuclear levels well known to occur for Ho 3+ cations, as 
discussed in Ref. [7^ and Ref. [73. This nuclear compo- 
nent was estimated by Blote et al. 76 for Ho2GaSb07. By 



subtracting it off from the (total) experimental specific 
heat, we can uncover the underlying magnetic contribu- 
tion and compare to the theoretically calculated Monte 
Carlo specific heat data, from which T pea k or C pcak can 
be determined directly (Fig.0. We note here, as recently 
observed in Ref.|71, that for Dy 2 Ti 2 07 there should be a 
hyperfine nuclear contribution to the specific heat man- 
ifesting itself at a temperature below T < 0.4 K in the 
data of Ref. if one uses the typical hyperfine contact 
interaction expected for a Dy 3+ insulating salt. The ab- 
sence of the high temperature I/T 2 tail of the nuclear 
specific heat (on the descending low temperature side of 
the magnetic specific heat of Dy2Ti207) below 0.4 K in 
Fig. 3a is, therefore, somewhat puzzling. 79 



□ Ho 2 Ti 2 G 7 experimental data 
— estimated nuclear component 
o magnetic (exp. - nuclear) 
MC magnetic component 
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FIG. 8: The total experimental specific heat of Ho2Ti2C>7 
is shown by the open squares. The expected nuclear contri- 
bution is indicated by the line, while the resulting magnetic 
specific heat estimation is shown by the open circles. Near 
0.7 K the estimation is prone to a large error. Dipolar spin 
ice simulation results are indicated by the filled circles. 

The shoulder-like feature in the estimated magnetic 
contribution to the experimental specific heat data of 
Fig. [3] (open circles) near 0.7 K can be entirely elim- 
inated by adjusting the nuclear hyperfine splitting by 
~ 2% percent around the value estimated by Blote for 
Ho2GaSb07, resulting in an exceedingly good agreement 
with the Monte Carlo results down to T = 0.4 K. Such 
a slight adjustment to account for any small deviations 
in the hyperfine parameters of 4/ rare-earth ions (depen- 
dent upon electric field gradients, chemical shift, etc.) 
would seem reasonable. However, we do not do this in 
order to emphasize that the unbiased use of the estimated 
nuclear specific heat contribution from the isostructural 
material Ho2GaSb07^ already allows for a very good 
agreement with the theoretical magnetic specific heat. 

Having determined J nn and D nn for Ho2Ti207 from 
specific heat measurements, we are able to compare the 
experimental clastic neutron scattering against that de- 
termined via the Monte Carlo simulations. The results, 
reported in Ref. |U show excellent agreement between 
experiment and simulation. More recent neutron scatter- 
ing experiments on the Ho2Sn207 show similar results^ 7 . 
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Such comparison between theory and experiments for 
Dy2Ti 2 07 is more difficult due to the large neutron ab- 
sorption cross section of naturally occurring Dy isotopes. 
Work in that direction using isotopically enriched sam- 
ples with 162 Dy isotope is in progress^ 

Numerical integration of the specific heat divided by 
temperature can be performed to determine the entropy 
of both Ho2Ti207 and Dy2Ti2C>7. Specifically, the en- 
tropy, S(T), removed between temperature T\ and T 2 , 
S , (T 2 ) — S(Ti), can be calculated using the thermody- 
namic relation: 

S(T 2 ) - g(Ti) = £ 2 ^ dT. (5) 

The results for Dy 2 Ti 2 7 arc illustrated in Fig.Eb- The 
entropy recovered between T = 0.4K, where C(T) is very 
small, up to a temperature T — 10K, is S(T = 10K) — 
S(T « 0) w 3.930 J mob 1 K _1 As we can see in Fig.Eb, 
the Monte Carlo data for S(T) at T = 10 K is slightly 
below the Pauling's value R{ln(2) - (1/2) ln(3/2)}. To 
perform the calculation of the recovered entropy between 
T = 10 K up to T — oo, we extrapolate the tem- 
perature specific heat C(T) for T > 10 K by match- 
ing the Monte Carlo value of C(T) at T = 10 K with 
the 1/T 2 high temperature paramagnetic temperature 
regime, C(T) = C^/T 2 for T > 10K. This gives a 
value Coo = 29.015 J mol -1 K, and an extra entropy 
of S(T = oo) - S(T = 10) = 0.145 J mol" 1 K _1 , hence 
a value S(T = oo) - S(T w 0) = 4.075 J mob 1 K" 1 , in 
exceedingly close agreement with Pauling's value, 4.077 
J mol^ 1 K _1 . Hence, we find that the simulation with 
the appropriate experimental parameters retains Paul- 
ing's entropy (Eq. ifTjl). similar to what is found experi- 
mentally for Dy 2 Ti 2 7 (Fig.Eb and in Ref.EiJ. A simi- 
lar experimental procedure was done using the magnetic 
contribution of the specific heat data of Bx>2Ti2 07, also 
giving a residual entropy close to Pauling's entropy^* 

While the above conventional Monte Carlo simulations 
of the model Hamiltonian for the spin ice compounds, 
Eq. (J2J), yields a reasonably successful quantitative the- 
ory of spin ice behavior in Ising pyrochlore materials, 
there still remains the second question (#2, Section ^Bj) 
as to why dipolar interactions, despite their anisotropic 
and long-range nature, do not (appear to) lift the macro- 
scopic degeneracy associated with the ice rules, and se- 
lect an ordered state. As a first attempt to address this, 
we investigate the spectrum of soft modes (i.e. critical 
modes or ordering wave vectors) accessible to the dipolar 
spin ice model within the context of mean-field theory. 



III. MEAN FIELD THEORY 

In this section we present the main results of a mean- 
field theory calculation aimed at determining the spec- 
trum of soft modes in the dipolar spin ice model. The 
details of the method can be found elsewhere 



The MC simulation results presented in the previ- 
ous section answer in the affirmative the question as to 
whether or not long-range dipole-dipole interactions in 
real materials are consistent with the manifestation of 
spin ice behavior in a temperature range < T < 9cw- 
However, these results do not address the question of 
whether or not a true ground state degeneracy and failure 
to order at any nonzero temperature is an exact symme- 
try consequence of the long-range dipolar interactions for 
(111) Ising spins on the pyrochlore lattice. A direct way 
to address this question is to ask whether or not there 
actually exists at the Gaussian level (i.e. MFT) a soft 
or critical mode in the model at a well-defined ordering 
wave vector q. The results of this calculation for finite 
distance cut-off, R c , of the dipole-dipole interactions have 
been reported in a conference proceedings We briefly 
review the essence of the calculation and extend it to 
untruncated (true long-range 1/-R 3 ) dipole-dipole inter- 
actions using the Ewald summation technique. In MFT 
the Ewald technique is implemented in q-space, in con- 
trast to real space for MC simulations. The approach is 
briefly discussed in Appendix lAl 

Our MF derivation begins with the Hamiltonian for 
(111) Ising spins, Eq. J5J, but expressed in terms of the 
Ising variables, af, and local quantization axes, z a , 

i,j a.b 

where 

J ah (i,j) = J(z a ■ z b )6 RaKR _ (7) 

(z«-z b 3(z°-R° fc )(z fc .R° b ) \ 
M [\R$\ 3 |R£T J' 

Recall that indices a and b denote the sub-lattice and 
R° b is the vector that connects spins af and er b . The 
pyrochlore lattice is a non-Bravais lattice which is de- 
scribed in a rhombohedral basis with four atoms per unit 
cell located at positions r a given by (0, 0, 0), (1/4, 1/4, 0), 
(1/4, 0, 1/4), and (0, 1/4, 1/4) in units of the conventional 
cubic unit cell of size a = i? nn v / 8- Each of these four 
points define a fee sub-lattice of cubic unit cell size a. 
Using H from Eq. ©, we form the free energy of our 
system, 

F = Tr{pH} + TTr{plnp}, (8) 

where p is the many-body density matrix. 

The first step in the mean-field approximation entails 
replacing p with the product of single-particle density 
matrices (p({uf}) = IL a Pi( a t ))• Next, a variational 
free energy is obtained by treating the p\ (af) as varia- 
tional parameters subject to the constraints Tr{pf} = 1 
and Tr{p"af } = m", where mj is the local magnetiza- 
tion, or order parameter. The resulting variational free 
energy is transformed to momentum space. Our conven- 
tion for the Fourier transform employs the position of the 



local magnetization, R", 



E'< 



(9) 



from which the spin interaction matrix can be expressed 
in terms of its Fourier components, 



iVcell „ 



(10) 



where iV cc n is the number of four atom unit cells, i.e., 
fee lattice points. We note that the above convention for 
the Fourier transform produces a symmetric 4x4 matrix 
i/(q) at all values of q. An alternate convention defining 
the Fourier transform with respect to the Bravais lattice 
points, Ri, results in a complex (Hermitian) J(q), as 
discussed in Ref. IH2I From Eqs. © and 1)10(1 . we write 
the quadratic part of the MF free energy, F^ 2 \ 

/ (2) (T)^EES{^ Qb -^ ab (q)}™-q . ( n ) 
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FIG. 9: X max (q)/D vs. q for q in the (00/) direction, in 
units of 2n/a, for various cut-off distances, R c , and for the 
infinite range limit simulated by the Ewald summation tech- 
nique. Calculations were made at J nn /D nn — on the spin 
ice side of the phase diagram, Fig. 



where f (2) (T) = F^(T)/N ccn and T is the temperature 
in units of 1/fce- Diagonalizing /' 2 '(T) requires trans- 
forming to the normal modes of the system, 

4 

m a = U a > a (q)<Z>« , (12) 

a=l 

where the Greek index a labels normal modes, {$q} are 
the amplitudes of the normal modes, and U(q) is the 
unitary matrix that diagonalizes t7(q) in the sub-lattice 
space, with eigenvalues A(q), 

t/t(q)J( q )[/( q ) = A(q) . (13) 

In component form, U a,a (q) represents the a-component 
of the a-eigenvector. We express (T) in terms of nor- 
mal modes as 

/ (2) (r) = i^^$«{r-A a ( q )}$« q . (14) 

q a 

In our approach, a minus sign was pulled out in front 
of the Hamiltonian in Eq. 6, therefore, an ordered state 
first occurs at the temperature defined by the global max- 
imum eigenvalue, 

T c = max q {A maa; (q)} , (15) 

where \ max (q) is the largest of the four eigenvalues (a = 
1, 2, 3, 4) at wave vector q, and maxq indicates the global 
maximum of the spectrum of X max (q) for all q. The 
value of q for which A"(q) is maximum is the ordering 
wavevectorciord ■ 

In Ref. 1^3) we evaluated J{q) in the spin ice regime 
of our model (J n n/D nn = 0) directly from the inverse 



transform of Eq. i|10[l for various real-space cut-off dis- 
tances of the dipolar term. For the rather dramatic ap- 
proximation of nearest-neighbor distances, we found a 
complete degeneracy (q-independence) of A max over the 
whole Brillouin zone. This corresponds to the q-space 
signature of the degenerate nearest-neighbor spin ice 
model of Harris and Bramwell, 21, and, equivalently, of the 
nearest-neighbor (global) Ising antiferromagnet model of 
AndersoniS As R c was increased, a i? c -dependent q or d 
appeared. In the limit of infinite range interactions, i.e., 
R c — > oo, we predicted q or d = (001) (in units of 2ir/a). 
Results for A max along the (00/) direction as a function 
of R c are shown in Fig. 

Treating dipolar interactions via the Ewald method, 
where the true long-range nature of the interactions are 
respected, we observe a completely smooth (without rip- 
ples) and quasi-degenerate soft mode spectrum with a 
global maximum (critical mode) at (001). The spectrum 
of A max (q)/D in the (hhl) plane of the pyrochlore lattice 
for J nn = is shown in Fig. 1101 while the results along the 
(00/) direction are included in Fig.^J Inspection of cor- 
responding eigenvectors of the doubly degenerate critical 
mode indicate a two-in two-out spin ice structure, where 
the spins on sub-lattices a = 1, 3 point opposite to those 
on sub-lattices a = 2,4 in a tetrahedral unit, see Table 
[I] Ref. |52j discusses how these soft modes can be used 
to reconstruct the long-ranged ordered (equal moment) 
structure, where the thermal local magnetization is the 
same on all sites, and which corresponds to the long- 
range ordered two-in two-out spin ice state of Fig. 4. 

With the implementation of the Ewald method, we are 
able to study the symmetry properties of the long-range 
dipolar interactions in a controlled manner. This is a de- 
sirable feature for a problem like spin ice because dipole- 
dipole interactions beyond nearest neighbor distances in- 
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FIG. 10: The scaled maximum eigenvalues, A maa; (q)/Z), in 
the (hhl) plane for the same spin ice model described in Fig. 03 
The dipole-dipole interactions are treated with the Ewald ap- 
proach. 



troduce perturbations onto the highly degenerate soft 
mode spectrum of the parent state, the nearest neigh- 
bor spin ice model. The form of the dipolar interaction, 
which couples spin and spatial degrees of freedom, per- 
mits both positive and negative contributions to the total 
dipolar energy as R c is increased, so the interactions are 
self-screening. However, this symmetry is not exact and 
the value for q or( j in MFT depends crucially but unpre- 
dictably on how far the sum is carried out. The variabil- 
ity in q or( j is readily observed for short cut-off distances, 
e.g., R c ~ lOOnn. However, subtle effects are still present 
even at very large cut-off distances and can produce an 
anomalous ordering wave vector. With the dipolar inter- 
actions cut-off at R c — lOOOnn, an incommensurate crit- 
ical mode is found at q ~ (0.025,0.025, 1.0). This mode 
is singly degenerate (unlike the doubly degenerate mode 
at (001) in the Ewald limit) and the eigenvectors do not 
predict a two-in two-out spin ice structure. We under- 
score, again, that this effect is subtle, as is demonstrated 
by the difference in the eigenvalues at R c = lOOOnn, 
A max (0.025, 0.025, 1.0) - A max (0, 0, 1) « 3 x 10~ 4 . 

In magnetic systems, the development of spin-spin cor- 
relations can be observed in the q-dependent suscepti- 
bility, x(q)- I n the Gaussian approximation one has 
X(q) oc (T - A"(q)) _1 (see Appendix 01 . Therefore, 
the critical mode that minimizes the free energy, /( 2 )(T), 
also controls the magnetic correlations as T — > T c , with 
T c given by Eq. (|15l) . In the paramagnetic (PM) regime, 
T » #cw, one expects all modes to contribute to x(q). 
This is also the regime in which MFT applies, T > T c , 
and is expected to provide quantitatively accurate re- 
sults. Therefore, a build up of PM correlations is under- 
stood in terms of an underlying critical mode that marks 
a transition to an ordered state at T c . This mean- field 



TABLE I: The two maximum eigenvalues and correspond- 
ing eigenvectors of ^(qord) a t q or d = (001) for spin ice with 
Ewald evaluation of the dipolar interactions. The U a ' a (q OT d) 
are normalized eigenvectors, where the weights indicate a two- 
in (positive value) two-out (negative value) spin ice structure. 
In particular, the spins on sub-lattices a — 1 and a = 3 point 
out of and spins on sub-lattices a — 2 and a = 4 point in to 
the tetrahedron. 

a A Q (q ord )/D tr- Q (q ord ) 

1 31575 ^(-1,1,0,0) 

2 3.1575 -W0, 0,-1,1) 



approach has been used to study the PM elastic neu- 
tron scattering in the (111) pyrochlores^For the dipolar 
spin ice model, one obtains results that are in agreement 
with the experiments and simulations of Ho2Ti2C>7 for 
R c > 250nm2£ The best results are found when dipoles 
are treated with the Ewald technique. If one employs too 
small a cut-off distance, e.g., much less than R c « lOOnn, 
then one finds that the PM scattering is concentrated in 
regions of the first zone that are inconsistent with experi- 
ments. This is especially true at cut-off distances studied 
m Refs.|3l] and|32j. Again, we believe this is a direct con- 
sequence of the failure of a finite dipolar sum to much 
restore the symmetry of the dipolar Hamiltonian, i.e., to 
produce a quasi-degenerate A max (q)-spectrum. 

From MFT, we can also determined the value of 
J nn /D nn at which the ordering changes from an "all- 
in— all-out" q = state (i.e., from large negative anti- 
ferromagnetic Jnn) to the (001) long-range ordered spin 
ice state. We find that the transition between the two 
states occurs at J n n/D nn = —0.905. This is in full agree- 
ment with the value found in Monte Carlo simulations 
results for the transition between all-in— all-out q = 
ordering and spin ice behavior4» 

Having obtained strong evidence from MFT that there 
exists a well-defined, unique ordering wave vector in 
the long-range dipolar spin ice model at the Gaussian 
level, we can proceed with our search for a dynamically- 
inhibited transition to long-range order in the model by 
the use of a conventional Metropolis single spin-flip MC 
simulation. 



IV. DYNAMICAL FREEZING AND LOOP 
MOVES IN MONTE CARLO SIMULATIONS 

A. Dynamical Freezing in Conventional Single Spin 
Flip Monte Carlo Simulations 

The mean field theory results presented in the previ- 
ous section make it clear that from the point of view of a 
strictly equilibrium (statistical mechanics) magnetic or- 
dering phenomenon, the dipolar spin ice model of Eq. (J2J) 
should be a rather conventional system with a unique and 
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well defined ordering wavevector and staggered magne- 
tization order parameter. The question then becomes: 
why don't Monte Carlo simulations of the dipolar spin ice 
model, or in fact the real spin ice materials themselves, 
develop the long-range ordered phase predicted by mean 
field theory? The problem turns out to lie in the local 
single spin flip dynamics employed within the Metropolis 
algorithm and, similarly, the local spin dynamics at play 
in the real materials. As we will show below, Monte 
Carlo simulations of the dipolar spin ice model using 
single spin flips experience a dynamical freezing at low 
temperatures. This arises due to the existence of large 
energy barriers separating distinct quasi-degenerate spin 
ice configurations, and prevents the simulation (and the 
real materials) from finding its true ener getic ally- favored 
long-range ordered ground state (see Ref.|82|for a related 
problem). 

Observation of the acceptance rate A(T) (percentage 
of accepted Monte Carlo steps) of the dipolar spin ice 
simulations makes it immediately apparent that out-of- 
equilibrium freezing occurs at low temperatures, that is 
below T ~ 0.4 K (as illustrated in Fig. Illh ) for the J nn 
and D nn parameters appropriate for Dy2Ti2 07. Fig. lllb 
shows that A(T) can be parametrized by a Vogel-Fulcher 
temperature dependence as found in numerous freezing 
phenomena: A(T) oc exp(A/(T — Tf reczc )), where the 
freezing temperature, Tf rcczc , is introduced in an ad hoc 
fashion. In Fig.lHb. A = 1 K. 

It is clear that in order to investigate the existence of a 
true energetically-favored ground state in the dipolar spin 
ice model, a standard Monte Carlo simulation employing 
local single spin flip dynamics is insufficient. Indeed, as 
Fig. El shows, these dynamical processes are frozen-out 
at T just slightly below 0.4 K for a model of Dy 2 Ti 2 7 . 
For J nn and D nn appropriate to describe Ho2Ti 2 07r^ 
the single spin flip Monte Carlo acceptance rate falls be- 
low 10 -6 at a temperature near 0.6 K. Without ascribing 
any deep significance to it, it is interesting to note that 
this freezing out in the simulation at 0.4 K and 0.6 K 
corresponds rather closely to the temperatures at which 
freezing is found in Dy2Ti207 41 and Ho2Ti207, 4 ° respec- 
tively. 

This freezing-out occurs due to large free energy bar- 
riers separating the (almost) degenerate ice-rules states, 
which develop rapidly at a temperature of order T pea k, 
and which are associated with introducing a single spin 
flip to a tetrahedron obeying the ice rules. As discussed 
above, and further supported by the mean-field calcu- 
lation, the effective (ferromagnetic) nearest neighbor in- 
teraction J c ff favors the ice-rules configuration. As the 
temperature drops the Boltzmann weight exp(— 4J c ff/T) 
suppresses the probability that a spin flip will take a 
given tetrahedron into an intermediate, thermally acti- 
vated non ice-rules obeying configuration. Thus single 
spin flip Monte Carlo moves are, for all practical pur- 
poses, frozen-out and dynamically eliminated within the 
simulation when T << J e ff. 
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FIG. 11: (a) Single spin flip Monte Carlo step acceptance 
rate A(T) for a simulation of Dy2Ti207. The simulation be- 
comes frozen when the acceptance rate falls to zero, (b) The 
logarithm of the acceptance rate plotted versus inverse tem- 
perature, minus some freezing temperature Tf ree ze, follows a 
Vogel-Fulcher type law. 



B. Loop Moves in Monte Carlo Simulations 

In order to explore the low temperature ordering prop- 
erties of dipolar spin ice, one needs a Monte Carlo 
algorithm with non-local updates that effectively by- 
pass the energy barriers that separate nearly degener- 
ate states and allows the simulation to explore the re- 
stricted ice-rules phase space that prevents ordering in 
the model, 58 59 In other words, we employ non-local dy- 
namical processes to restore ergodicity in the Monte 
Carlo simulation, and then use this new algorithm to ex- 
plore and characterize the long-range ordered state that 
arises out of ice-rules manifold and which is energetically 
favored by the dipolar spin ice model. 

We first identify the true zero energy modes that can 
take the near-neighbor spin ice model from one ice state 
to another exactly energetically degenerate ice state. An 
example of these zero modes, or loops, is shown in Fig. 1. 
We take as an initial working hypothesis that in the 
dipolar spin ice model, with interactions beyond near- 
est neighbor, the system freezes into an ice-rules obey- 
ing state. This is indeed what we found: in all of the 
tests we performed, systems simulated using conventional 
single spin-flips always froze-out in an ice-rules obeying 
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state with no "defects" (by defects we mean violations 
of the Bernal- Fowler ice rules). With interactions be- 
yond nearest-neighbor, these loop moves become quasi- 
zero modes that can take the dipolar spin ice model from 
one ice-rules state to another without introducing spin 
defects into tetrahedra in the lattice. This allows all of 
the quasi-degenerate spin ice states to be sampled er- 
godically, and facilitates the development of a long-range 
ordered state by the system at low temperatures. 

Within the Monte Carlo simulation, we use the 
Barkema and Newman^i loop algorithm originally de- 
signed for two dimensional square ice models, and adapt 
it to work in a similar manner on the three dimensional 
pyrochlore lattice. In the context of square ice, we tested 
two types of loop algorithms, the so-called long and short 
loop moves. In the square ice model, each vertex on 
a square lattice has four spins associated with it. The 
vertices are analogous to tetrahedron centers in the py- 
rochlore lattice. The ice rules correspond to "two spins 
pointing in, two spins pointing out" at each vertex. In 
the Newman and Barkema algorithm, a loop is formed 
by tracing a path through ice-rules vertices, alternating 
between spins pointing into and spins pointing out of 
the vertices. A "long loop" is completed when the path 
traced by the loop closes upon the same spin from which 
it started. A "short loop" is formed whenever the path 
traced by the loop encounters any other vertex (tetrahe- 
dron) already included in the loop - excluding the dan- 
gling tail of spins (Fig. 1121) . 




FIG. 12: Long and short loops formed by the Newmann and 
Barkema algorithm—*— on a square ice lattice. Vertices are 
represented by points where lattice lines cross. Each vertex 
has two spins pointing in and two spins pointing out, however 
for clarity, only spins which are included in the loops are 
shown. Starting vertices are indicated by large black dots. 
On the left is an example of a long loop, which is completed 
when it encounters its own starting vertex. On the right is 
a short loop, which is complete when it crosses itself at any 
point. Dark gray lines outline completed loops. The excluded 
tail of the short loop is shown in light gray. 

We now generalize the Barkema and Newmann loop al- 
gorithm for our study of the three dimensional pyrochlore 
lattice spin ice problem. In this system, the smallest com- 



plete loop that is a zero mode on the pyrochlore lattice 
consists of six spins (see Fig. QJ. Such a loop was pre- 
viously identified by Bramwell and Harris^i and also by 
Anderson^ in the context of the spinel lattice. However 
using the above loop algorithm, much larger loops are 
possible. When used with the pyrochlore lattice (Fig.QJ, 
such a loop must pass through two spins on each tetrahe- 
dron. A loop always "enters" a tetrahedron through an 
inward pointing spin, and "leaves" a tetrahedron through 
an outward pointing spin. The periodic boundary con- 
ditions of the lattice may also be traversed with no ill 
consequences. If we form a closed loop in this manner, 
and each spin is reversed on it, the entire system stays 
in an ice-rules state. However, small dipole-dipole en- 
ergy gains or losses may be procured due to small energy 
differences between the old and the new ice-rules state. 
These small energy changes caused by the loop moves are 
evaluated via a Metropolis algorithm within the Monte 
Carlo. Specifically, a loop move that takes the system 
from one ice-rules state to another one of lower energy 
is automatically accepted, while a loop move that takes 
the system to a higher energy ice-rule state (with energy 
difference 5E between the two states) is accepted with 
exp[— 5E/{k^T)\ > rnd, where rnd is a random number 
taken from a uniform distribution between and 

Before use in a full-scale Monte Carlo simulation, the 
long and short loop algorithms are subjected to a vari- 
ety of characterizing tests on the three dimensional py- 
rochlore latticed The first test is a study of the rela- 
tive speed (measured by CPU time) of the al gori thms 
for different sized lattices. As reported in Ref. 159. it is 
found that the small loop algorithm creates loops that 
approach a finite size limit as the system size increases. 
The long loop algorithm continues creating larger and 
larger loops, that scale approximately linearly with the 
number of spins in the simulation cell. This forces the 
algorithm to become drastically slower for the larger sys- 
tem sizes considered. 

Second, tests are carried out to investigate how the two 
different loop algorithms handle defects that break the ice 
rules on a tetrahedron. As we know, above the "spin ice 
peak" in the specific heat of the dipolar spin ice model, 
the ice rules (two spins pointing into a tetrahedron, two 
spin pointing out) are generally not obeyed. However, to 
retain detailed balance, we want our loop algorithm to at- 
tempt to form loops at temperatures above the onset of 
spin ice correlations. The attempt to create a loop is sim- 
ply aborted in the case where the loop path encounters a 
defect (either a three-in one-out vertex, or an all-in or all- 
out vertex) . The simulation does not attempt to flip any 
spins on an aborted loop, and the Metropolis al gori thm 
is not employed in this case. As reported in Ref. l59l the 
result of ice rules defects on loop algorithm performance 
is significant. Within the long loop algorithm, the in- 
clusion of only one defect spin per one thousand spins 
in the system causes almost half of all loops which are 
attempted to be aborted on the grounds that they have 
encountered the defect, with efficiency decreasing drasti- 
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cally as more defects are included. In contrast, the short 
loop algorithm remains 87% efficient with the inclusion 
of one defect in one thousand, and retains an efficiency 
that is much better than the long loop algorithm as more 
defects are present in the systemiS 

We use both algorithms to perform a true finite tem- 
perature Metropolis Monte Carlo simulation of the dipo- 
lar spin ice model. We code both algorithms separately, 
and run simulations using the regular procedure (cooled 
slowly from a high temperature, equilibrating carefully at 
every temperature step). Our preliminary Monte Carlo 
results for both the short and long loop are given in 
Fig. Even though the data in this figure has low 
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FIG. 13: Preliminary data for the low temperature magnetic 
specific heat (a) and energy (b) of the dipolar spin ice Monte 
Carlo, system size L=3, with simulation parameters set for 
Dy2Ti2 07. The data represent an average taken over ap- 
proximately 10 5 production Monte Carlo steps. Closed circles 
are data from a simulation of the short loop algorithm, open 
squares are data obtained using the long loop algorithm. Low 
temperature features are discussed in the next section. 



to choose between the two based solely on their perfor- 
mance properties measured above. As alluded to above 
and detailed in Ref. |5lJ the short loop algorithm works 
more efficiently within the requirements of our simula- 
tion. However, the disadvantage with using the short 
loop algorithm is that each loop does not cover as large 
of a percentage of spins within the lattice. It is not clear 
to us, without investigating the computational perfor- 
mance (i.e. autocorrelation times) of both algorithms in 
much more quantitative detail, whether a small number 
of long loops is better at bringing the system to equilib- 
rium than a larger number of short loops for a fixed CPU 
time. However, with the additional observation that the 
long loop algorithm can pass over itself a number of times 
during its creation, effectively losing additional efficiency 
in this manner, we ultimately choose to perform the ma- 
jority of the simulations on the dipolar spin ice model 
using the short loop algorithm. 

With the short loop algorithm chosen for the simula- 
tions, we re-investigate the Monte Carlo Metropolis al- 
gorithm acceptance rate. Since each loop successfully 
created (i.e., not aborted by encountering an ice defect) 
by the short loop algorithm is still subject to rejection 
by the Metropolis condition on the basis of its change in 
system energy, we expect the maximum acceptance rate 
to be somewhat less than the maximum efficiency of the 
algorithm given above. Results for the loop acceptance 
rate are shown in Fig. [21 Clearly, the loop algorithm be- 
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FIG. 14: The Monte Carlo acceptance rate for the short loop 
algorithm (circles) as compared to the single spin flip algo- 
rithm (line), for a simulation using Dy2Ti2 07 parameters. 
Acceptance rates were calculated as percentages of attempted 
Monte Carlo steps (MCS). 



statistics (only 10 5 Monte Carlo production steps per 
spin), it is clear that both the short and long loops pro- 
mote roughly the same thermodynamic behavior in the 
Monte Carlo simulation. The low temperature features 
(specific heat peak and energy discontinuity at T w 0.2 
K) of Fig.^Jare induced in the same manner by both al- 
gorithms. These features will be discussed in much more 
detail in Section V. Since both the long and short loops 
display equivalent results in the Monte Carlo, we are free 



comes effective in the temperature range where the single 
spin flip algorithm looses its ability to explore all possible 
configurations of the system. Above about 1 K, the num- 
ber of accepted loops is very low, due to the fact that the 
system is not entirely in an ice-rules configuration. As 
the simulation is slowly cooled, ice-rules constraints be- 
gin to develop, and the loop algorithm begins to work 
efficiently, moving the system between different ice-rules 
states. In Fig. ^[ a sharp drop is observed in the loop 
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acceptance rate at approximately 0.18 K. As discussed 
in the next section, this corresponds to the temperature 
where a phase transition develops in the system, which 
locks the system into a long-range ordered state and elim- 
inates Monte Carlo dynamics once and for all. 



V. LOOP MONTE CARLO INVESTIGATION OF 
THE TRANSITION TO LONG RANGE ORDER 

As suggested by the results above, the short loop algo- 
rithm is successful in restoring ergodicity in the simula- 
tion. As a consequence of this, we observe a low temper- 
ature phase transition in the model. The most familiar 
and robust indicator of a thermodynamic phase transi- 
tion (as opposed to dynamical freezing) is a finite size 
remnant of a singularity (divergence or discontinuity) in 
the specific heat at the transition temperature T c . In our 
simulation of Dy2Ti207, a sharp cusp in the specific heat 
is observed at a temperature below the spin ice peak (See 
Fig. 115b ). The feature in the specific heat and the abrupt 
drop in energy at the same temperature gives good pre- 
liminary evidence that the loop algorithm is successful 
in allowing a phase transition to occur at a temperature 
of T c = 0.18 K. This is the same temperature where the 
loop algorithm acceptance rate goes to zero in Fig. ITU 
The energy curve shows a discontinuous drop at T c (e.g., 
latent heat) for large lattice sizes, suggesting a first order 
phase transition. In the remainder of this section, we at- 
tempt to characterize this ordered state, and the phase 
transition that leads to it. The first step is to identify 
the order parameter associated with the low temperature 
ordered state. 

Direct inspection of the spin directions at T < 0.18 
K reveals that the ordered state is a long-range ice-rules 
obeying state with zero magnetic moment per unit cell 
and commensurate with the pyrochlore cubic unit cell 
(see Fig. . This state corresponds to the critical mode 
found above in the mean-field calculation. There are 
twelve symmetrically equivalent spin configurations for 
the ground state as explained below, two for each cu- 
bic axis direction and their spin reversed states. The 
ordering wavevector q lies parallel to one of the cubic 
axis directions, specifically q = (0, 0, 2n/a) or one of its 
symmetrically equivalent directions. To construct the or- 
dered state, first consider a starting tetrahedron with 
its six possible ice-rules states. For a given ordering 
wavevector q, this tetrahedron selects one of the four 
possible spin configurations (two independent configura- 
tions and their spin-reversals, Sf — > — S") with a total 
magnetic moment for the tetrahedron perpendicular to q. 
The entire ordered state may then be described by planes 
(perpendicular to q) of such tetrahedra. The wavelength 
defined by this q physically corresponds to antiferromag- 
netically stacked planes of tetrahedra, which means that 
a given plane has tetrahedra of reverse spin configuration 
to the plane above and below it. 
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FIG. 15: The low temperature magnetic specific heat (a) and 
energy (b) of the dipolar spin ice Monte Carlo, system size 
L=4, with simulation parameters set for Dy2Ti207. Closed 
circles are simulation data run with the short loop algorithm, 
open triangles are data obtained using the single spin flip 
Metropolis algorithm. In the inset of (b), the energy shows 
an apparent discontinuity at a critical temperature T c ~ 0.18 
K. The broad feature in the specific heat at T » 1 K indicates 
the rapid development of the spin ice rule obeying states. The 
sharp feature at T c is the appearance of a phase transition 
to a ground state being made (dynamically) accessible via 
the non-local loop dynamics. Note that these results are of 
higher statistics than those for Fig. 1131 specifically, 1 x 10 5 
equilibriation and and 3 x 10 J production Monte Carlo Steps 
were used. In addition, the system size is increased to L = 4 
as opposed to L = 3. Also, note that the location of the 
specific heat peak is at roughly the same temperature and is 
narrower than for L = 3, indicating a finite-size effect on the 
singular behavior of C(T). 



We construct the multi-component order parameter 



\b m = 

a N 



N/4 4 
j=l a=l 



(16) 



This type of labeling is natural given that the pyrochlore 
lattice can be viewed as an FCC lattice with a "down- 
ward" tetrahedral basis (Fig. [TJ . Thus j labels the FCC 
lattice points of the pyrochlore lattice, and the index a 
sums over the four spins comprising the basis connected 
to each j. The index a labels the three possible symmetry 
related q ordering wavevectors. For a given q a , as de- 
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scribed above, there are two ice-rules configurations and 
their reversals which can each form a ground state. Thus 
rn = 1.2 labels these possibilities with the phase factors 
{4>™}, describing the given configurations m. Each Ising 
variable ct° has a value +1 or -1 when a spin points into 
or out of its downward tetrahedron j, respectively. 

As written in Eq. (|16fl . St™ has six degenerative com- 
ponents, each of which can take on a value between 
and 1. Upon cooling through the transition, the system 
selects a unique ordered configuration, causing the corre- 
sponding component of '3/™ to rise to unity and the other 
five to fall to zero (provided the finite size system is sim- 
ulated over a time scale less than the ergodic time scale 
where full spin symmetry is restored). The component 
which rises to unity is equally likely to be any one of the 
six, selected by random through spontaneous symmetry 
breaking. 

Fig- El is a plot of (\&) for three system sizes, where 



w = 



\ 



E E (*s 

m=l a=l 



(17) 



is the magnitude of the multi-component order param- 
eter. These three curves illustrate important finite size 
effects for (fy). For T < T c the different lattice sizes 
produce identical order parameters. By contrast, (\&) 
for the smaller lattice size displays pronounced rounding 
near T c and an increased residual value for large T . The 
larger lattice size produces an order parameter with a 
clear discontinuity at T c . This discontinuity in the order 
paramater combined with the discontinuity of the total 
energy in Fig. 15b can be viewed as strong preliminary 
evidence for a first order transition. 




FIG. 16: The q = (0,0, 2n/a) order parameter. Curves are 
shown for system size L = 2, L = 3, and L= 4. 



We now argue the need to study this phase transition 
with greater numerical accuracy. This is necessary par- 
tially to confirm rigorously its first-order nature. More 
importantly, once this is done, we want to use the data 
to confirm the full recovery of Paulings entropy through 
an estimation of the latent heat released by the tran- 
sition. To begin, we note that there are a number of 
criteria at one's disposal to demonstrate the occurence 



of a first-order transition in a Monte Carlo simulation. 
In particular: 

1. The order parameter (^S>) should have a clear dis- 
continuity at T c . 

2. The energy probability histogram, H(E), should 
have a double peak at T c , which identifies the co- 
existence of two distinct phases at T c . 

3. There should be a latent heat at the transition, 
identifiable by a discontinuity in the internal energy 
for large system sizes. 

4. In the Monte Carlo, the height (maximum) of the 
specific heat, C pca k, and the magnetic susceptibil- 
ity, Xpoak, should be proportional to the simulation 
volume: 



Cpeak i Xpeak OC a + bL d 



(18) 



where a and b are constants and d is the dimension 
of the lattice, in our case equal to three (d = 3). 

5. The minimum of the fourth order energy 
cumulantpS 



V = l 



should vary as 

^min =V + cL- d 

where V ^ 2/3. 



(19) 



(20) 



6. The temperature T peak (L) at which C pea k or x pea k 
have a maximum should vary with the simulation 
volume as: 

T peak (L) - T c + cL- 1 (21) 
where c is a constant, and T c = T pcak (L — > oo). 

We have already confirmed the first condition of our 
list. To check for the second condition the energy proba- 
bility histogram was calculated by binning the simulation 
energy values for every Monte Carlo step as the system 
passes through the transition from higher to lower tem- 
peratures (Fig. IT7I). Above T c , we observe a single peak 
Gaussian-like distribution of energies. At T c , the energy 
probability distribution shows a double peak, character- 
istic of the coexistence of two phases found at a first 
order phase transition. Below T c we would normally ex- 
pect to see a Gaussian peak. However, in our case, the 
histograms below T c are distorted by the accumulation of 
energies into the lowest bins due to the proximity of the 
transition to the ground state. At zero temperature, we 
expect all of the energies to lie in the bin for the lowest 
energy. 

The next condition on our list is the observation of a 
latent heat at the transition. Fig. El shows the energy 
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FIG. 17: The energy probability distribution histogram for 
three temperatures: T = 0.178 K, T = 0.180 K (T c ), and T = 
0.182 K. For T > T c (filled circles), a single peaked Gaussian 
histogram is present. At the transition temperature (hashed 
rectangles), a second peak appears which has a lower mean 
energy. As the temperature falls below T c (filled triangles), 
the peak with the higher mean energy disappears, and the 
system energy eventually gathers in the lowest bin. 
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FIG. 18: Details of the simulation energy near the transition 
for different system sizes. 



near the transition for three different system sizes. A 
clear discontinuity develops as we increase the system 
size. The energy discontinuity, AE, as read off of this 
graph for L = 4, is 

AE ps 0.248 J mol -1 . (22) 

This behavior is also consistent with the transition being 
first order. Below, we use this AE value in calculating 
the entropy recovered at the transition (entropy jump). 
Another calculation of the latent heat at the transition 
comes from the finite size scaling of condition 4 abovesS 

where L d is the system volume as before, and a is the 
intercept of the graph of C pea k(£) vs. L to test for con- 
dition 4, which we now discuss. 



When attempting to quantify the relationships in con- 
ditions 4, 5 and 6 on our list, we notice a problem. The 
extremely sharp nature of the transition makes accurate 
estimates for these quantities almost impossible using a 
traditional temperature cooled MC simulation of the hy- 
brid single spin flip— loop algorithm. The reason is that 
the transition temperature region is so narrow, and first- 
order metastability effects are so strong, obtaining accu- 
rate data for quantities such as Cp Ca k or V m i n very near 
to T c is extremely difficult. As shown in Fig. El the en- 
ergy probability histogram near a first order transition 
displays a double hump. The energies that occur be- 
tween these humps correspond to system configurations 
that are strongly suppressed by the Boltzmann probabil- 
ity distribution near the transition. We call these "in- 
terface configurations" i24 Traditional Monte Carlo sim- 
ulations try to "avoid" these interface configurations as 
the system is cooled through the transition, because of 
their suppression by the Boltzmann factor which is the 
basis of the Metropolis condition. Therefore, the simula- 
tion often behaves poorly in this region, moving quickly 
through interface configurations to find more favorable 
configurations nearby in configuration space. This can 
lead to erratic behavior and poor statistics in thermo- 
dynamic quantities of interest near the phase transition, 
thereby reducing the numerical accuracy of the quantities 
used in finite-size scaling. 

To overcome this problem, Berg and Neuhaus^i pro- 
posed the multicanonical method, which is designed to 
enhance configurations that have energies which occur 
between the double hump of the probability distribu- 
tion. If these interface configurations are artificially en- 
hanced, the simulation does not avoid this energy range 
as strongly and better statistics can be obtained. The 
version of the multicanonical Monte Carlo algorithm that 
we use is that proposed by Hansmann and Okamoto, 85 
originally developed to be used in the context of protein 
folding simulations. The core of the method is: 

Perform Monte Carlo simulations in a multicanoni- 
cal ensemble instead of the usual canonical ensemble. 
Then, obtain the relevant canonical distribution by us- 
ing the histogram reweighting techniques of Ferrenberg 
and Swendsenm^ From this, calculate the thermodynamic 
quantities of interest. 

In the multicanonical ensemble, we define the proba- 
bility distribution by 

Pmu (E) = 9(E)w um {E) = congtant (24) 

where g{E) is the density of states, w mn (E) is the mul- 
ticanonical weight factor (not temperature dependent), 
and Z mu is the associated partition function. The distri- 
bution is constant, meaning that all energies have equal 
weight, which sometimes leads to the name "flat his- 
togram" method. This flatness is important because it 
ensures that configurations in the interface region of the 
transition are not suppressed. 

Unlike for the canonical ensemble, the multicanonical 
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weight factor w mn is not a priori known. This turns out 
to be the crucial step of this scheme: finding an accurate 
estimator of w mu that makes the distribution p mu (E) flat 
over the energy range of interest. The details of how to 
do this are somewhat involved, and will not be explic- 
itly outlined here. The reader is referred to the relevant 
technical references for detailsi£2i2& Our procedure fol- 
lows that of Ref. |8^ very closely. 

Assuming that we can find a good estimator of w mu , 
our method proceeds as follows: 

1. We find an accurate estimator of the multicanoni- 
cal weight factor so that p mu (E) is reasonably flat 
over an energy range that includes the transition 
interface. 

2. With this weight factor we perform a multicanoni- 
cal simulation at one given temperature T slightly 
higher than T c . 

3. During this simulation run, we gather statistics for 
the physical variables of choice (for example, the 
energy E). These variables are weighted according 
to the multicanonical distribution. 

4. From this single simulation, we then obtain the 
Boltzmann-distributed variables at any tempera- 
ture for a wide range of temperatures using a 
reweighting technique. 

We use the reweighting technique proposed by Fer- 
renberg and SwendsenjSk which allows us to transform, 
or reweight, data obtained from another distribution (in 
our case the multicanonical distribution) to the relevant 
Boltzmann distribution, at some inverse temperature (3. 
We use this to obtain an estimate for a given physical 
quantity in the canonical distribution. 

We collect data within the multicanonical distribution 
and use it to calculate the specific heat for the dipolar 
spin ice model. For the smallest system size considered, 
we accurately reproduce the specific heat over the tran- 
sition using the flat histogram method. Fig. I19fa shows a 
comparison between the specific heat of an L = 2 system 
obtained using the histogram method at one tempera- 
ture, and the traditional Monte Carlo procedure with 
8 x 10 5 equilibriation steps and 2 x 10 6 data production 
steps for every temperature point. The CPU time that 
it took to get the histogram data was a small fraction 
of the time it took to obtain the regular Monte Carlo 
data. Fig. 119b is a similar result for the next highest sys- 
tem size, L = 3. The traditional Monte Carlo data was 
taken with 5 x 10 5 equilibriation steps and 1 x 10 6 data 
production steps. The histogram data was obtained in 
the same amount of time as for the L = 2 data, and it 
was only slightly more difficult to find a good estimate 
for w mu (E). The poor quality of the traditional Monte 
Carlo Metropolis data for L = 3 stands in stark contrast 
to the smooth data obtained using the multicanonical 
simulation. 
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FIG. 19: Specific heat curves over the transition temperature, 
for L = 2 (a) and L = 3 (b) system sizes. Closed circles repre- 
sent data obtained with Ferrenberg and Swendsen's histogram 
reweighting technique. Open triangles represent data taken 
using a traditional temperature cooled Monte Carlo simula- 
tion. 



Unfortunately, one difficulty with the multicanonical 
algorithm used here is that, in general, as the system 
size is increased, it becomes increasingly difficult to ob- 
tain a good estimate for a w mn (E) that would give a flat 
p mn (E). The critical temperature, T c , of the transition 
seems to be the quantity most sensitive to variations in 
the flatness of p mu (E). In contrast, the height of specific 
heat peak is fairly accurately determined for simulation 
sizes L = 2, 3, 4, and 5 (Fig. 12(1 [I . showing only a weak 
sensitivity to the flatness of p mu (E). 

The aforementioned error associated with T c for the 
L = 4 peak, as determined from simulations, is of the or- 
der of 0.04 K, and becomes increasingly more drastic for 
the larger system sizes. The variation in the height of the 
specific heat was found to be much less. Nevertheless, to 
combat any minor variation in peak height and obtain 
an accurate finite size scaling results, a statistical aver- 
age was done on several (~ 10) multicanonical weighting 
factors to obtain values for C poa k- These results are plot- 
ted in Fig. |2J A straight line fit to the data using linear 
regression gives 

Cp 0ak = 0.8924L 3 - 3.149. (25) 

The L 3 dependence of C pea k(£) shows that the finite size 
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FIG. 20: Specific heat of the transition to long-range order, for 
system sizes L = 2, 3, and 4. It was found that peak heights 
for these system sizes did not vary much with the flatness of 
Pmu(B). Critical temperatures T c were more sensitive to the 
flatness of the multicanonical distribution, and hence were 
harder to estimate than C pea k- 
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FIG. 21: Finite size scaling fit for the specific peak heights of 
the ordering transition. Data points represent the mean C pea k 
value for a given L. Error bars show one standard deviation. 



scaling is consistent with that expected for a first order 
transition. Also, as a second estimator of the latent heat, 
we can use the slope of this line and Eq. H23[) to extract 
AE. Doing so we deduce a latent heat of 

AE = 0.245 J mol -1 , (26) 

consistent to within 1% with the value obtained in Fig. 1181 
(for the L = 4 system) from reading directly off of the 
energy graph (see Eq. I^jll. 

This completes our study of the nature of the ordering 
transition in dipolar spin ice. As we have shown, the dis- 
continuity in the order parameter, the release of latent 
heat, the double peaked energy probability distribution, 
and the finite size scaling of the specific heat peak all give 
consistent and compelling evidence for the transition be- 
ing first order. As the technical details concerning this 
transition are understood, we can proceed to study where 
it, and the long-range ordered state which results from 



it, stand in our broader picture of ground state entropy 
found in experiments and in standard single spin flip sim- 
ulations of the dipolar spin ice model. Since we have 
confirmed the first order nature of the transition, the 
configuration of the ordered state, calculated the latent 
heat, and have reliable data for the specific heat through 
the transition, we are in a position to re-calculate the to- 
tal entropy that the dipolar spin ice model releases as it 
is cooled to low temperatures. This calculation must be 
done carefully. We know that in an infinite system, a first 
order transition is characterized by a cusp in the specific 
heat. If the transition is temperature driven, as in our 
case, this first order singularity is the latent heat. For 
an infinite system going through a first order transition, 
thermodynamics gives 

AS = F" ^dT+ r^>dT+—, (27) 
J T J T + T T ' y ' 

where AE/T is the latent heat contribution to the en- 
tropy (see Fig.UBfl. and T~ and and T+ are the temper- 
ature limits asymptotically close to T c , below and above 
T c , respectively (see Fig. ITS)) . 

To estimate a value for the entropy, we consider the 
system size L = 4 which has good statistical data for the 
widest temperature range. We integrate the low temper- 
ature data for the specific heat in Fig. [201 divided by tem- 
perature obtained from the histogram reweighting tech- 
nique (up to T « 0.21). For T > 0.21 K we use our regu- 
lar temperature cooled Monte Carlo data (canonical loop 
+ single spin flip) for the integration above this point, 
and up to 10 K, giving S(T = 10) - S(T « 0) = 5.530J 
mol -1 K _1 . To integrate up to T = oo, we follow 
the same high temperature extrapolation procedure de- 
scribed in Section II, giving S(T = oo) - S*(T = 10) = 
0.145J mol -1 K _1 . Doing this simple calculation, we find 
a total recovered entropy of 

S(T = oo) - S(T rs 0) = 5.675 J mol -1 K _1 , (28) 

less than 2% below our expected value of i?ln2 = 5.764. 
The inset of Fig. [221 clearly shows the entropy recovered 
by the low temperature transition. 

By considering the entropy recovered by the integra- 
tion of the finite size system specific heat over the transi- 
tion (inset, Fig. I22[). we confirm that it is approximately 
equal to the value of the entropy we would expect to re- 
cover from the latent heat of an infinite system. Using our 
latent heat calculations above (Eqs. 11' I'D and l|26|l '). this 
value is approximately AS Tc = AE/T C w 0.2465/0.180 = 
1.37 J/mol K, in good agreement with the jump in S(T) 
in the inset of Fig. |2U 

Taking all indicators together, we have demonstrated 
here that the transition to long range order at 180 mK 
recovers all residual Pauling entropy of the dipolar spin 
ice model. Thus we can assert that the degeneracy as- 
sociated with the spin ice model, and the corresponding 
value of zero point entropy, is lifted due to perturba- 
tions beyond nearest neighbor dipole-dipole interactions, 
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FIG. 22: The entropy calculated from integrating the simu- 
lated specific heat, as explained in the text. The entire value 
of entropy for the system (i?ln2) is recovered in the high 
temperature limit. The inset shows the details of the entropy 
recovered by the transition to long-range order. 



if equilibrium can be maintained at sufficiently low tem- 
peratures. 

We comment here on a detail we neglected to discuss 
above regarding the relationship between the mean-field 
theory and the results from the loop Monte Carlo simula- 
tions. In the Gaussian mean-field theory presented above 
the calculation that was performed was in effect an iden- 
tification of the soft-mode against which the paramag- 
netic phase becomes locally unstable upon cooling. The 
Monte Carlo simulation finds, however, that the thermo- 
dynamic transition to that ordered state is actually first 
order and occurs before the supercooled critical temper- 
ature is reached. 

To summarize our results for this section, we refer the 
reader to the dipolar spin ice Monte Carlo phase dia- 
gram, Fig. EI As illustrated there, the transition between 
the spin ice phase (which retains Pauling's entropy) and 
the q = (0,0, 2ir/a) ordered phase is independent of the 
strength of J nn . This is consistent with our understand- 
ing that the long range order results from perturbative 
interactions beyond nearest neighbor, caused by the long- 
range dipolar interaction. This is also what mean-field 
theory finds in the spin ice regime (J nn /Am > —0.905). 
We find that this first order line also slightly runs up the 
boundary between the antiferromagnetic ordered phase 
and the higher temperature paramagnetic phase, and 
that a tricritical point separates these two regions of the 
line, and occurs near the value J n n/D nn ~ —1.1. 

Due to the near-vertical nature of the phase bound- 
aries in this region, simulations run at a finite T and 
varying J nn help better map out the low temperature 
phase lines of interest. However, using this method, we 
observed that the simulations could easily get "stuck" 
in the previous spin configuration (either spin ice dis- 
ordered, q = (0,0, 2n/a) or AF q = 0) when cross- 
ing the vertical phase boundary. This history depen- 
dence is illustrated in the phase diagram as hystere- 



sis at low temperatures, mainly between the long-range 
ordered q = (0, 0, 2n/a) and antiferromagnetic q = 
phases. Regardless of this difficulty, we have confirmed 
from direct Ewald energy calculations at zero tempera- 
ture that the true zero-temperature phase boundary be- 
tween the q = (0, 0, 2n/a) and the AF q = phases lies 
at J nn / D nn — —0.905, in agreement with the result found 
above in the mean field calculations. 



VI. DIPOLAR SPIN ICE IN MAGNETIC FIELD 

A very interesting problem that pertains to dipo- 
lar spin ice materials is their behavior in an ex- 
ternal magnetic field, h. A number of recent 
experimentsSSiSiiSiSMii have shown a rich variety of new 
behavior when spin ice materials are subjected to such a 
field, which warrants some theoretical investigationifiiS 
Although not all of the relevant experiments can be de- 
scribed in this short section, we briefly outline some of 
the most important, referring the reader to the bibliog- 
raphy for further details on their methods and results. 

The first experiments on spin ice materials in an ap- 
plied magnetic field were performed by Harris et al.w& 
In a neutron scattering experiment, they applied a mag- 
netic field of strength 2 T along the [110] direction of a 
single crystal of Ho2Ti207, and looked for signs of or- 
dering. They found scattering intensity features which 
suggest evidence of two ordered magnetic structures, the 
so-called q = and q = X phases (Fig. |23J). As we 
will see below, these ordered structures are of fundamen- 
tal importance in our study of the ground states of the 
dipolar spin ice model. 

A quite interesting set of experiments was performed 
by Ramirez et (dM- and Higashinaka et al.fi& who sub- 
jected polycrystalline samples of Dy2Ti207 to a variety 
of different field strengths. Ramirez presented evidence 
of field-dependent phase transitions in a powder sample, 
manifested as sharp features in the specific heat at 

1. 0.34 K for h > IT 

2. 0.47 K for IT ~ h ~ 3T 

3. 1.12 K for h ^ 

where we have used h to represent the magnitude of 
the applied field h. Higashinaka and co-workers repro- 
duced the basic features of Ramirez' results down to 
T = 0.38 K, 88 confirming the existence of the two higher- 
temperature peaks only. The search for a microscopic 
explanation of these three peaks has been a driving force 
behind much of the experimental and theoretical work 
in this field over the past few years, and we discuss it 
further in the work that follows below. 

Another significant experimental study is the measure- 
ment of the single crystal magnetization curves (M vs. 
h) for the spin ice materials. Fukazawa et al. per- 
formed a number of experiments on single crystals of 
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FIG. 23: The pyrochlore unit cell projected down the z-axis. 
The symbols + and - represent the z component of the spin 
"head". Configurations are (a) q = and (b) q = X. The 
large arrows point in the (a) [100] and (b) [110] directions, 
and are included to aid in the discussion of ground states in 
the next section. Note that q = X, although similar, is a 
spin ice state distinct from the q = (0, 0, 2n/a) ordered state 
shown in Fig. 01 Specifically, the chains of spins parallel to 
the [100] direction are staggered antiferromagnetically in the 
zero field q = (0,0, 2n/a) ground state (Fig. 4), while they 
are ferromagnetically correlated and parallel to the field in 
the q = X state of (b) above. 



Dy"2Ti207, obtaining magnetization curves for the differ- 
ent applied field directions and a range of temperatures^ 
They showed that magnetization data at 2 K was consis- 
tent with the behavior predicted by the spin ice model, in 
particular the limiting field (large h) values of the mag- 
netic "anisotropy" (which we illustrate below). Very re- 
cently, measurements^*^ of the magnetization curves for 
h//[lll] (read "h parallel to the [111] crystal direction") 
have uncovered a novel macroscopically degenerate state 
corresponding to ice-like behavior on the kagome planes 
in the pyrochlore lattice*^! 

We take into account the applied magnetic field h in 
the dipolar spin ice model with a simple term added to 
the Hamiltonian (Eq. ©), 



ff' = -5>.s 4 a = -]T(h. 



(29) 



We work strictly with a classical Ising model and neglect 
any transverse field effects and perturbative changes to 
the moments arising in a strong field. In this classical ap- 
proximation, the field h couples to the spins through the 
simple scalar product Eq. I|29|) . That is, we neglect small 
corrections to the energy coming from the very small, 



though finite, local susceptibility perpendicular to the 
(111) direction. As well, we neglect quantum mechanical 
transverse- field effect that would arise from admixing the 
doublet ground state wavefunctions with that of the ex- 
cited crystal field levels. For Ho2Ti20y and Dy2Ti2C>7, 
the first excitation gap is (very roughly) A ~ 300 K. 28 
For the magnetic moments of approximately 10 [Ab for 
both Ho 3+ and Dy 3+ , this means a ground state Zee- 
man energy splitting of 12.8 K/Tesla. One can therefore 
safely neglect magnetic field, exchange and dipole-dipolc 
induced admixing for fields less than 10 Tcslas, assuming 
the worse case scenario where the excited doublet also 
split by about 10 K/Tesla. 

To gain a theoretical understanding of the experimen- 
tal behavior mentioned above, several insightful calcula- 
tions are possible, using only this simple classical Hamil- 
tonian and a knowledge of the possible ground states of 
Fig. [23] First, a geometrical understanding of how the 
magnetic field couples to classical spins on the pyrochlore 
lattice is desirable. We expect that application of a mag- 
netic field along the three principle symmetry axes of the 
crystal will result in different spin-field coupling behav- 
ior. To explore this, we begin by considering the non- 



interacting limit [h — > oo or J nn , D n 



0). In this 



case, the only constraints on the spins is the local (111) 
anisotropy and the coupling with the magnetic field. We 
can gain more insight by viewing a projection of a tetra- 
hedron down the cubic z-axis as in Fig. [^J 
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FIG. 24: A single tetrahedron projected down the 2-axis. 
Field directions are (a) h//[100], (b) h//[110], (c) h//[lll], 
depicted by the large arrow outline. Small arrows represent 
dipole moments coupled to the field. Empty circles represent 
decoupled spins. 

For h/ /[100], all four spins on a given tetrahedron are 
coupled with the field (i.e. all four have a non-zero dot 
product in Eq. (|29[l ). The expected lowest energy config- 
uration in the absence of spin-spin interaction is the one 
where all spins have their [100] components aligned with 
the field. Knowing this, we can calculate the h — > oo 
value for the M vs. h curves by considering the average 
moment M (in units of Bohr magneton per rare earth 
ion, ^s/R 3+ ), in the direction of the field h//[100]: 
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(30) 



Also, to calculate M in the appropriate units (/is/R 3+ , 
measured in experiments) one must include the factor fi, 
which is the magnetic moment of the appropriate rare 
earth ion. If R is Dy 3+ or Ho 3+ , fi w lOfiB- Note that, 
from Fig. [21k . this lowest energy spin-field coupled state 
is compatible with the ice rules. If we decorate the entire 
lattice with tetrahedra such as this, we recover the q = 
state of Fig. 123k . This suggests that this ordered state 
should be one of the ground states for the interacting 
dipolar spin ice model, with a sufficiently strong exter- 
nal field h//[100]. Indeed, this order has been observed 
experimentally^ on samples of Dy2Ti2C>7. 

For h//[110], only two of the four spins on a tetrahe- 
dron couple to the field. One expects that, with precise 
alignment of the sample, these other two spins would re- 
main decoupled even in the application of high magnetic 
fields. These decoupled spins are thus free to choose an 
ordering pattern that satisfies their dipolar interaction. 
Because of the complexity of the dipolar interaction, the 
ground state spin configuration is not immediately obvi- 
ous from studying the geometry. However, one expects 
any zero-temperature phase to be consistent with the ice- 
rules (see Fig. 124b). In the limit of very high applied field 
and perfect sample alignment, one expects the magneti- 
zation to approach 



M (h -> oo) 
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fj, S 0.4082 • fi. 
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Finally, for h//[lll], all four spins on a tetrahedron are 
coupled to the field. An interesting complication arises 
in this case due to crystal geometry; any high-field phase 
of the material will be inconsistent with the ice-rules, 
and the spins will form a three-in one-out (or its spin 
reverse) tetrahedral configuration (Fig. 124b). For zero 
temperature, both the long-range ordered ice-rules state 
and the three-in one-out state will exist for different field 
strengths. For low magnitudes of h, we expect a compe- 
tition between the exchange, dipolar and magnetic field 
parts of the Hamiltonian. At low enough temperatures, 
one predicts£2i£iiSi§2i£! that a plateau will develop in the 
magnetization curve due to the tendency of each tetra- 
hedron to stay in the ice rules up to a critical field. If we 
couple three of the spins to the magnetic field, and leave 
one to oppose the field but obey the ice rules, we find a 
magnetization of 

M (h= "small") 

= 171 + [Tll] + [lTl] + [TTl]) ' 7! ' M 



— • fi = 0.3333 • /i. 



(32) 



above, the last spin vector [TTl]) to break the ice rules, in 
favor of minimizing its energy with respect to the field. 
In this case the high field magnetization is 
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In the limit of very high applied field, we expect the 
spin that is coupled anti-parallel to the field (in the case 



We find that our Monte Carlo is successful in repro- 
ducing the high field limiting values of the experimental 
M vs. h curves*^ In addition, we find that the Monte 
Carlo also reproduces the plateau expected for h//[lll] 
and low temperatures However, because these large h 
results are easily obtainable for a nearest-neighbor spin 
ice model (D nn — * 0), we won't discuss them further in 
this work. The reader is referred to Ref. |5^ and Ref. |(5(] 
for the detailed results of this study. 

A numerical calculation of interest that is easily per- 
formed is the Ewald energies of the various ground state 
configurations that we have encountered so far in the 
dipolar spin ice model. The spin ice configurations that 
we consider are both the q = and q = X phases iden- 
tified by Harris, 20 and the q = (0,0, 2ir/a) ground state 
identified previously in this work (Fig. In addition, 
we expect a "three-in one-out" state to become the low- 
est energy state for some critical field along the h//[lll] 
direction. Figs. 1251 and 1261 are the results of these ground 
state energy calculations for a system size L=2 and 
parameters appropriate for Ho2Ti207 (J nn = — 0.52K, 
D nn = 2.35K). As expected, we find the same qualita- 
tive behavior for calculations involving Dy2Ti20y param- 
eters. In addition, because all of the ground state config- 
urations considered are commensurate with the unit cell 
of the pyrochlore lattice, these calculations scale trivially 
for sizes L > 2. 

Fig. 12 5b confirms that the q = configuration becomes 
the lowest energy state for large field strength (h > 0.034 
T) for h//[100], as expected from simple geometrical con- 
siderations. Recall that in the h//[110] case, there exist 
two decoupled spins per tetrahedron, and subsequently 
no lowest-energy configuration is obvious from the geo- 
metric field coupling Eq. (|29(l . However, one may antic- 
ipate that for h//[110], the decoupled spins (Fig. 121b 1 
order in "chains" perpendicular to the [110] direction, 
arranged in such a way as to partially satisfy the dipo- 
lar interaction. We find that this is precisely the q = X 
state, which Fig. 125b shows to be the lowest energy state 
at for h > 0.023 T. 

Consequently, Fig. 125b has direct relevance to exper- 
iments by Fennel et al. and Hiroi et al. that were per- 
formed on Dy 2 Ti 2 7 with a magnetic field h//[110]^a 
Fennel et al. observed neutron diffraction patterns that 
showed Bragg scattering at q = "points" and diffuse 
scattering at q = X "points" , but no true q = X long- 
range order. They suggested that this behavior would 
arise from long-range ferromagnetic order occurring in 
field-coupled spin chains (called a chains by Hiroi§£), 
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FIG. 25: The T=Q energies per spin of the three ice-rules 
ordered states of the dipolar spin ice model, as a function 
of applied internal field h — |h| along the (a) [100] and (b) 
[110] directions for J nn = -0.52 K and D nn = 2.35 K (i.e. 
Ho2Ti2C>7 parameters). 



and short-range "antiferromagnetic" order occurring in 
the field-decoupled spin chains (/? chains)^ In this ar- 
gument, the true ground state is a q = X structure that 
is dynamically inhibited from being accessed on experi- 
mental timescales. 

Specific heat measurements by Hiroi et al. were used 
to extract the specific heat contributions of both the a 
chains and the (3 chainSfSS, They suggest that the spe- 
cific heat due to the f3 chains resembles that which one 
would expect for a low-dimensional spin system without 
long-range order. They also argue for the presence of 
geometrical frustration in the triangular sub-lattice that 
contains the j3 chains. If such a frustration exists, it 
might be expected to destabilize the "antiferromagnetic" 
correlations between these chains that would otherwise 
lead to q = X order. Therefore, Hiroi et al. argued 
against true long-range order for the system, rather that 
the chains become effectively isolated and behave as 
"pure" one-dimensional ferromagnetic systems without 
long-range order in the ground state. 

At this point, Fig. 125b is consistent with the idea of 
long-range q = X order for the dipolar spin ice model 
with a magnetic field h//[110]. As we will discuss be- 
low, finite-temperature Monte Carlo calculations on the 
dipolar spin ice model also support the idea that, similar 
to the development of q = (0, 0, 27r/a) order in the zero- 
field case, the development of q = X order for h//[110] 
may in some cases be dynamically inhibited in experi- 
mental systems with local spin dynamics. The failure of 
the frustration of the j3 chain sub-lattice invoked by Hi- 



roi et al. 62 to destroy the long-range q = X order may 
be another example of the small energy scale left over 
by the infinite-range dipole-dipole energy (Eq. J2J), and 
why such interactions must be handled carefully using 
techniques such as the Ewald method. 

It should also be noted here that the q = and q = X 
lines are parallel in Fig. 125b only for samples that are 
perfectly aligned with h along the [110] crystal axis. This 
is an important phenomenon one must consider when 
comparing theory and experiment, as only a small crystal 
misalignment will partially couple spins on the (3 chains 
to the field. Because precise alignment of a crystal is 
often very difficult, the possibility of misalignment of the 
order of a degree must be taken into consideration when 
studying single crystal data with h//[110]. Repeating 
our ground state energy calculation for misalignment of 
one degree along the [100] direction, one finds a crossing 
of the q = X and q = lines in Fig. 125b at about 1.3 
T, the q = configuration being of lowest energy above 
this field strength^ 
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FIG. 26: The T=0 energies per spin of the three ice-rules or- 
dered states and the three- in one-out spin state of the dipolar 
spin ice model as a function of applied internal field h along 
the [111] direction. The q = X state becomes the ground 
state at 0.029 T. The three-in one-out configuration becomes 
the ground state at around 1.4 T, breaking the ice rules for 
each tetrahedron. The q = X line is always 0.534 J mol -1 
below the q = line at any given field. 



Also, as Fig. [SB] confirms, the three-in one-out spin 
configuration become the lowest energy state for large 
h//[lll]. Interestingly, the q = X state is the ground 
state for 0.029 < h < 1.4T. This is consistent with the 
idea that that the z component of the field for the [111] 
direction gives a zero net Zeeman contribution to the 
unit cell for both the q = and q = X spin con- 
figurations. Therefore, the only energy scale difference 
left over between the two states comes from the [110] 
component of the field, meaning that q = X will be 
slightly lower in energy than q = (as in Fig. 25b) for 
all h//[lll]. Most other features of the h//[lll] field 
(such as the intermediate-field plateau and the high-field 
breaking of the ice-rules) are readily explainable in a 
nearest-neighbor spin ice model without dipolar interac- 
tion, hence we will discuss them no further in this work. 
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We now turn to a set of preliminary results of finite- 
temperature Monte Carlo simulations for the dipolar spin 
ice model in a magnetic field. As mentioned previously, 
the desire to explain the three field-independent specific 
heat peak a 29 i 88 of polycrystalline Dy2Ti20y has driven 
much of the experimental and theoretical interest in this 
field over the past few years. Ramirez et al. were the first 
to suggest^ 9 - that some of these peaks can be attributed to 
a fraction of crystallites whose [f 10] axes happen to align 
(closely) with the applied magnetic field. More generally, 
we can interpret this argument as saying that magnetic 
moments which are not strongly coupled to the magnetic 
field through Eq. l(2"5|) are free to contribute to a dipolc- 
induced phase transition (and therefore sharp peaks in 
the specific heat) at low temperatures. Historically, the 
h//[110] coupling of Fig. 124b was considered to be the 
most likely scenario to provide these field-decoupled spins 
in a finite number of crystallites in the polycrystalline 
sample. 29 However, one may in fact argue that crystal- 
lites with only one field-decoupled spin would occur in 
much greater number in a real polycrystalline sample. 
This is due to the fact that, for a given crystallite orien- 
tation, there are an infinite number of applied magnetic 
field directions for which a given sub-lattice is decoupled 
(one of these being h/ /[112] 6 — ), corresponding to a rota- 
tion degree of freedom in the choice of h that does not 
exist in the two spin field-decoupled case (h//[lf0]). 

Hence, we carry out finite-temperature Monte Carlo 
simulations on the dipolar spin ice model for various field 
directions to look for signs of an ordering transition in 
the specific heat. In the case of h//[100] and [110], the 
ground state is known and hence an order parameter can 
be constructed, facilitating the identification of any pos- 
sible phase transitions. We perform simulations on the 
spin ice model with the added term Eq. I|29|l . employ- 
ing both single spin flips and loop moves (and ignoring 
the demagnetization effects discussed in Appendix C). 
For fields parallel to the [110] crystal axis, of magnitude 
large enough to favor the q = X ground state, but still 
relatively small, we found that the simulations were able 
to find this fully ordered state only when the loop moves 
were employed (see Fig. 1271 . This can be understood by 
studying the structure of the q = X state and its coupling 
to this field direction, as shown in Figs. 123b and!24b. As 
discussed above, the field-decoupled spins occur in long 
chains (f3 chains) in the bulk material. Within a single 
[3 chain, spins tend to point ferromagnetically along one 
direction as dictated by the exchange and magnetic field 
energies. In addition, the dipolar term weakly couples 
nearby chains, and in order to minimize this coupling 
energy, neighboring chains can be expected to seek out 
a unique ordering pattern. However, to explore different 
chain-orientations at lower temperatures, whole chains 
must be flipped at once in the simulation (imagine go- 
ing from the q = state to the q = X state in Fig. 1231 . 
Thus, if the system is attempting to settle into the q = X 
state at a temperature well below the onset of ice-rules 
correlations, these global chain-flip moves must be em- 




0.0 0.5 1.0 1.5 2.0 

T (K) 



FIG. 27: The specific heat of Dy 2 Ti 2 07 for h//[f 10], obtained 
using simulations with single spin flips and loop moves. Inset: 
the q = X order parameter calculated by the simulations. 
The single spin flips (SSF) are unable to find the true ground 
state of the system, while the loops moves (Loops) allow the 
system to order into the q = X structure. The large value of 
the high-temperature tail of the order parameter is a finte-size 
effect. 



ployed in the simulation. Luckily, such chain-flip moves 
are a subset of the generic loops created in the Monte 
Carlo loop algorithm (due to periodic boundary condi- 
tions), and thus no major modification of the simulation 
procedure is needed. 

As we see in Fig. there is evidence of a feature 
in the specific heat that corresponds to the temperature 
where the q = X order parameter jumps to essentially 
the saturation value of one. The feature in the specific 
heat and the corresponding jump in the q = X order 
parameter are at approximately 0.3 K, show that this 
is not the same transition as the q = (0,0, l-nja) tran- 
sition of the previous section. For small applied field, 
the transition temperature depends on the strength of 
the applied magnetic field. For example, at 0.075 T the 
transition temperature moves up to 0.4 K, and requires 
the loop moves to be manifest. For fields of the order of 
0.01 T and larger, the transition has risen to high enough 
temperatures that single spin flip dynamics are still suf- 
ficiently present to promote development of the ground 
state, without help from the loop moves. Details of this 
are illustrated in Fig. |^H1 Note that the phase transition 
illustrated in these two figures appears to be strongly 
first order, and hence obtaining good error control is dif- 
ficult with the current statistics (8x 10 4 equilibration and 
8 x 10 4 production MC steps). A more detailed study em- 
ploying multicanonical techniques is currently in progress 
to look more closely at the exact temperature and field 
dependence of this new phase transition. 

This sharp specific heat feature persists to a magnetic 
field value of h ~ 0.15 T. At higher field strengths, the 
singularity broadens and gets absorbed into the spin ice 
peak of the specific heat at T ~ 1.2 K. However, for all 
h > 0.5 T, the order parameter is observed to jump to 
unity at around 0.9 K, which appears to be the saturation 
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FIG. 28: The specific heat of Dy 2 Ti 2 7 for h//[110], obtained 
using simulations with single spin flips and loop moves. Inset: 
the q = X order parameter calculated by the simulations. 
Both the single spin flips and the loop moves are able to find 
the true q = X ground state of the system. 



value for very strong fields. Hence, although the nature of 
the transition to long-range order changes at some field 
h ~ 0.15 T, the fact remains that the simulations can 
always find the q = X ground state up to infinite field 
values. 

For larger field values (h > 0.2T) and h//[110], we re- 
produce the experimental- 2 , features of the specific heat 
for Dy2Ti207 nicely. Fig. [55] which shows our results 
for < T < 8 K for several field values, is relegated to 
Appendix [U] where we discuss demagnetization effects. 
Comparing this figure to Fig. 2, Ref. [62, shows that our 
model quantitatively reproduces the development of two 
peaks as h is increased above 0.5 Tesla. The T ~ 7 
K peak, clearly visible for h = 1.0 T (Fig. H^l), is a 
Schottky-type peak which corresponds to the freezing of 
field-coupled (a) spin chains that are along [110]. 62 The 
lower peak at T ~ 1 K is attributed to the development 
of long-range correlations between /^-chains, however its 
exact relationship to the ordering of the model as ob- 
served by the jump in the q = X order parameter is still 
under study. 

Finally, we performed finite-temperature Monte Carlo 
simulations on the dipolar spin ice model for fields paral- 
lel to the [100], [111] and [112] crystal directions. In par- 
ticular, we used field strengths larger than those which 
would favor q = (0,0,27r/a) order (see Figs. 1251 and I26[) . 
In the case of h/ /[100], we see the gradual development 
of q = long-range order as the simulation cools down 
through the spin ice peak, and no sharp singularities 
in the specific heat which may signal a phase transi- 
tion. Simulations for h//[lll] also shows no unexpected 
anomalies in the specific heat. Interestingly, work on the 
h//[112] case has not yet produced any signs of an or- 
dering transition in the field-decoupled spin sub-lattice, 
as one might expect from our discussion above. Work 
addressing this problem is ongoing. 

To summarize this section, we have performed several 



calculations of the properties of the dipolar spin ice model 
in an external magnetic field. Ewald energy calculations 
of various spin configurations reveals the preferred T = 
ordering for various field directions. In particular, for 
h//[100] the q = structure is the ground state, while 
for h//[110] we find the q = X state becomes energet- 
ically favored. For finite-temperature Monte Carlo sim- 
ulations of the dipolar spin ice model, we find that the 
system settles into these ground states for the respective 
field directions, although non-local dynamics are needed 
to find the q = X state for h//[110] at low fields. We 
also predict a sharp first-order phase transition to the 
q = X state and a corresponding specific heat spike for 
0.02 <h< 0.15 T in Dy 2 Ti 2 7 . Interestingly, unlike the 
transition to the q = (0,0,27r/a) state in zero field, the 
transition to the q = X state can be found with single 
spin flips in the Monte Carlo for fields 0.10 < h < 0.15 
T, although it is unknown whether local dynamics will 
be manifest in the corresponding range in a real experi- 
ment. Finally, for h > 0.15 T and parallel to [110], the 
simulation exhibits the broad peaks in the specific heat 
observed in experiments^ 6 ^ It is likely that the T = 1.1 
K specific heat feature corresponds to the T = 1.1 K fea- 
ture found in polycrystalline samples of Dy2Ti207i22iS 



VII. CONCLUSION 

We have reviewed much of the early experimental and 
theoretical work on the static magnetic properties of spin 
ice. We have also clarified our principle point of view that 
long-range dipolar interactions are consistent with and 
responsible for the physics observed in spin ice materials 
based on Dy 3+ and Ho 3+ rare earth ions. Support for 
our perspective resides in the detailed Monte Carlo and 
mean-field calculations presented in this paper. 

Monte Carlo simulations were performed on the dipo- 
lar spin ice model with the long-range dipole-dipole in- 
teractions treated via the Ewald method. Using a single 
spin flip Monte Carlo method, we were able to study the 
development of the spin ice manifold. We found that 
spins freeze out at temperatures O(IK) with a macro- 
scopic degeneracy or residual Pauling entropy. We also 
found that single spin flip dynamics are not effective at 
equilibrating the system, thus making it impossible to 
determine the ordered state of spin ice by this technique. 

Mean-field theory (Gaussian approximation) was ap- 
plied to the same dipolar spin ice Hamiltonian with the 
dipolar interactions treated via the Ewald method in q- 
space. There, we showed that an ordering wave vec- 
tor may be selected and that a proper treatment of the 
long-range dipoles is crucial to achieving a consistent 
picture with the experiments. A key point is that the 
symmetry of self-screening is not exact for the dipolar 
Hamiltonian. In the end, we found a quasi-degenerate 
spectrum emerges with a commensurate critical mode 
(q = (0, 0, 2tt/o)) and a two-in two-out spin ice struc- 
ture. 
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In order to find the ordered state of spin ice in a 
Monte Carlo simulation, we developed a non-local algo- 
rithm that employs loop moves (or updates) when in- 
side the spin ice manifold. These loop moves represent 
the "nearly" zero energy collective dynamics that allow 
our model to sample the highly degenerate phase space 
of spin ice. Application of this method at temperatures 
within the spin ice manifold, i.e., T ~ IK, leads to the se- 
lection of a single spin ice ground state configuration with 
q = (0,0, 2tv/o). The loop Monte Carlo and mean-field 
results agree. In addition, we find a first order transition 
to the ground state at T c w 180 mK, which recovers all 
of the residual Pauling entropy of the spin ice manifold. 
Our physical understanding of spin ice is aided by the 
picture that any collective dynamics in real spin ice ma- 
terials are inhibited by a freezing process as the system 
enters the temperature range where the ice-rules fulfilling 
manifold develops, i.e., Tfrcezc ~ 0.4 K for Dy2Ti207 and 
Tfroozo « 0.6 K for Ho 2 Ti 2 7 compared to T c « 180 mK. 

On the strength of the experimental evidence and the 
success of the dipolar spin ice model, we assert that both 
Ho2Ti207 and Dy2Ti20y are spin ice materials. 

Finally, we have reflected on the application of a mag- 
netic field to the spin ice materials as means of exploring 
the possible structures of the spin ice manifold and to 
further characterize the interactions present in these in- 
triguing systems. We find excellent agreement between 
the dipolar spin ice model and many experimental stud- 
ies to date. In addition, we find evidence for a low- 
temperature ordering transition to a q = X ground state 
for small magnetic fields parallel to the [110] crystal di- 
rection, that has not yet been observed. Some of the 
results presented here regarding the behavior of spin ice 
are intriguing. This argues for more theoretical, numeri- 
cal and experimental work, to resolve all the perplexing 
issues at stake. 
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APPENDIX A: EWALD 

We give only a brief overview of the Ewald^ technique 
as it applies to dipole-dipole interactions in Monte Carlo 
simulations and at the mean-field level. A more detailed 



discussion of the method can be found in Refs. [73j and 
I90L The mean-field case as it applies to moments on the 
pyrochlore lattice is treated in depth in Ref. |H3 

The dipole-dipole interaction is an infinite sum that 
falls off as the inverse cube of the separation distance 



between dipoles, 1/R! 



afc|3 



Hence it is a conditionally 



convergent series. The point of the Ewald method is to 
convert this slowly converging lattice sum into of two 
absolutely (rapidly) converging series, one in real space 
and the other in Fourier space. The general lattice sum 
for (111) Ising dipoles on the pyrochlore lattice is 



A = 




3(z a -R%)(z b -R%) 
|R?*I 5 



i,j a, 6 



tj ab 



(Al) 



x=0 



where the spin variables af have been dropped for nota- 
tional convenience. The dipole sum excludes terms with 



R 



ab 



0. Absolute convergence is forced on the sum 



inside the curly brackets of Eq. (|A1|) by use of a conver- 
gence factor. The form of this convergence factor differs 
depending on whether the dipolar sum is performed on 
N particles in real space (e.g., Monte Carlo and molecu- 
lar dynamic simulations) or in the thermodynamic limit 
in momentum space (mean- field theory). 

In our work, MC simulations are performed on 3 di- 
mensional lattices of LxLxL cubic cells of the pyrochlore 
lattice under periodic boundary conditions, thus there 
are N — 16 xLxLxL spins in the simulation cell. The 
separation of moments within a simulation cell is given 
by Rfj . The dipolar energy for any pair- wise interac- 
tion is calculated within the minimum image convention 
by summing replicas of the iV-site simulation cell over 
spherical shells of radii n = L(n x ,ri v ,n z ) (n x ,n y ,n z are 
integers) with the inclusion of a spherical convergence 
factor e -s|n| . The effect of the convergence factor is 
removed from the final form of the Ewald equations by 
imposing the limit s — > 0. Therefore, the starting point 
for the Ewald method is the dipole-dipole pair interac- 
tion, 



E 



,-s|n| 2 



n + R° 



x=0 

(A2) 



where J2 n r means that n = is omitted whenever R?- 
0. The point charge distribution, l/|n + R? fc — x|, 
rewritten with the aid of the T-function identities, 

= -4 [°°t-i/ 2 e-W 2 dt 



Using Eq 





3j , the pair- wise interaction becomes 
1 



-(z a -V x )(z b -V x )- 



dt 



(A3) 
(A4) 

(A5) 
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and the remainder of the Ewald for calculation for Af b 
follows arguments outlined in Ref. |9(J. The Ewald equa- 
tions for a Monte Carlo simulation can also be found in 
Appendix A of Ref. |6^ Af b is calculated for each pair- 
wise interaction, {{i, a), (J, 6)}, in the simulation cell, but 
this need be done only once because the spins are fixed to 
the lattice points. These pair interactions are stored in 
a look-up table and used in the stochastic sampling and 
measurement procedures of a Monte Carlo simulation. 

In mean-field theory, one considers the Fourier trans- 
form of the dipole-dipole interaction; therefore, the term 
e -iq-R i:j - g mc i uc } ec [ m Eq ||Aljl and plays the role of a 
convergence factor as discussed in Ref. With a con- 
vergence factor that is periodic as opposed to spherical in 
R, the derivation of the q-dependent Ewald equations fol- 
lows the method of long waves introduced by Born and 
Huang, Ref. As noted in Section IIIII the momen- 
tum dependent dipole-dipole interaction is determined 
between sub- lattice sites and is a component of j7 oh (q), 
where J7(q) is 4 x 4 matrix. Using Eq. (|A4|) . we can write 
the dipolar contribution to l / afc (q) as 

A a '»(q) =-(z a - V x )(z 6 -V x )^= / dt (A6) 

V 71 " Jo 

Jx=0 

where means the i = j term is excluded when 

a = b. The remaining steps in the derivation of the q- 
dependent Ewald equations can be found in Ref. I52L In 
this work, „4° ,& (q) was determined in at each q-point in 
the first zone of the (hhl) plane of the pyrochlore lattice 
and used in the formation of J7(q) in Eq. 1)11(1 . We note 
that the value of .4 a ' b (q) in the limit q — > dependents 
on direction. The value of A a ' b (0) can be related to the 
demagnetization factor i 91 ! 92 



APPENDIX B: x(q) FROM A HIGH 
TEMPERATURE SERIES EXPANSION 

We demonstrate that the mean field description of the 
static susceptibility, 46 x(q), can also be derived from a 
high temperature series expansion (HTSE) to lowest or- 
der in p. In our HTSE, where there is no reliance on a 
mean-field approximation, we show that the same func- 
tional dependence on the eigenvalues that leads to a min- 
imum in the free energy also defines the maximum in the 
response function of the local magnetization, which we 
interpret as a transition to a long-range ordered state 
at the Gaussian level. We use the model for local (111) 
Ising spins on the pyrochlore lattice, Eq. but the final 
results are applicable, with minor modifications for the 



spin components, to a general Heisenberg Hamiltonian. 
The Ising (111) Hamiltonian is 

H = -\HH^ ab ^» ( B1 ) 

i,j a,b 

where erf = ±1 are Ising variables and J ah (i,j) is given 
by Eq. 0. 

The q-dependent susceptibility is defined as the 
Fourier transform of the two-point correlation function, 

= ^EE^" • * b )W ■ <uV q - R?J \ (B2) 

-^*cell , . . 

where (3=1 /T with temperature in units of ks ■ We note 
that at this point x(q) is a 4x4 matrix for Ising spins on 
the pyrochlore lattice. A HTSE for (erf cr 6 ) is expressed 
most clearly as an expansion of cummulants^ 

°° (_a\m 

<°M> - E -Or w*J Hm ^ ( B3 ) 

m— 

On the left hand side of Eq. i|B3|) . (...) represents a ther- 
mal average and is thus a trace over states at finite T. 
The cummulants, (...) CJ are expressed as traces over the 
T = states (i.e., (...) = Tr{...}/2 N ). For our pur- 
poses, we need consider only the first two terms in the 
expansion, 

(ofo*) « «^) c - 0{o?o*-H) c . 

In the case of Ising spins, a non-zero trace has an even 
number of spin variables at each site. The zeroth order 
cumulant is determined trivially, 

(<7>*> c = ( O f *) = 5^0. 

The first order contribution involves two terms, 

{ay)H) c = (oy)H) - «cr}) (tf) , 

but the second does not contribute because (H) a cx 
( a i a m)o = 0. Therefore, one has 

«er^) = -^^J rf (l, m )«^>f)» 

Ira cd 

= -J a \hj)- 

The q-dependent susceptibility in the high temperature 
limit reads, 

M J- E £(*° • ^)(M a " + PJ ab (hj))e^. 

^*ccll , - ■ 
a,b 1,3 

Performing the Fourier transform of J ab {i 1 j), the inverse 
of Eq. (JTUJ, we obtain 

X (q) = [3j2(Z a ■ z h W h + ft7° 6 (q)), (B5) 

a . b 
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where J{— q) = J(c\) for a symmetric interaction ma- 
trix. The spin-spin interaction matrix and thus the sus- 
ceptibility are diagonalized via the normal mode transfor- 
mation given by Eq. The unitary matrix, U(q), that 
diagonalizes ^7(q) contains the orthonormalized eigenvec- 
tors of J(q), i.e., J2 a C/ a ' Q (q)C/ b ' Q (-q) = I, where I is 
the 4x4 unit matrix. We use this to rewrite the two 
terms in Eq. I|B5|) . 



(B6) 



* af, = ££t/ a ' a (q)!7 b ' a (-q) 



a a,b 



and 



J ab {d) = A Q (q)[/ a ' Q (q)[/ b < Q (-q) 



(B7) 



The expression for x(q) now becomes, 
X (q) = I>° ' f + /3A Q (q))t/ a ' a (q)C/ fo ' Q (-q). 



a a.b 



(B8) 

In the high temperature limit, (3 — > 0, and 1 + /3A Q (q) « 
1/(1 — /3A"(q)); therefore, the static, q-dependent sus- 
ceptibility now reads, 



x(q) = /?EE 



Q a. 6 



(z Q • z b )[/ a < a (q)[/ b > Q (-q) 
(l-/3A°(q)) 



(B9) 



= a sr IE.^ n,a (q)l 2 

(l-/3A-(q)) • 

Hence, as T approaches the ordering temperature defined 
by Eq. lfT5)l. T c = A max (q ord ), the susceptibility diverges 
and signals a transition to a long-range ordered state. 



APPENDIX C: DEMAGNETIZATION EFFECTS 

When doing finite temperature Monte Carlo simula- 
tions on magnetic materials in an applied magnetic field, 
the effect of the boundary of the simulation cell must be 
carefully considered. For systems of interest, the dipolar 
spin ice Hamiltonian is augmented with a field depen- 
dent term, Eq. I|29|l . The inclusion of this term in our 
Monte Carlo simulations leads to subtle effects. In a mi- 
croscopic Hamiltonian, the field h referred to in Eq. H29|) 
is the sample internal field, i.e., the magnetic field that 
directly couples to each magnetic dipole moment. How- 
ever in real materials, bulk demagnetization effects alter 
the magnitude of the internal field in a complicated man- 
ner that depends on sample size, shape, alignment, and 
surrounding medium. In general, experimentalists define 
three separate quantities (the magnetic flux density B, 
the magnetic field strength H, and the magnetization M) 
to account for these effects. In a macroscopic material, 
these quantities are related by 



where B is the independent quantity controlled in the 
experiment, but H is the field strength that couples to the 
spins (through which the bulk susceptibility is defined). 
In order to benchmark an experiment to a theory such 
as ours, an attempt must be made to relate the external 
experimental (applied) B controlled by external sources 
of current to the internal h of our Hamiltonian. From an 
experimental side, this amounts to knowing the internal 
M associated with the specific sample being measured. 
This M is in general not easily deduced; although, for 
certain sample shapes (e.g., ellipsoids of revolution) it 
is at least uniform, and fairly accurate estimates can be 
made. The procedure of correcting for M to obtain H is 
called making a demagnetization correction. 

Theoretically, demagnetization effects are incorporated 
into a Monte Carlo simulation by imposing certain 
boundary conditions on the microscopic Hamiltonian in 
question. As described earlier, we use the Ewald sum- 
mation method to calculate the long-range dipolar inter- 
actions of our model. We follow the standard approach 
in which the pair wise interactions are evaluated by sum- 
ming over periodic copies of the N site simulation cell 
until convergence is obtained^ effectively simulating the 
infinite range nature of the dipoles. A consequence of this 
technique is that the finite size nature of the simulation 
cell is suppressed. We are, therefore, faced with the ques- 
tion of how to interpret an "infinite boundary" . If one 
wishes to simulate materials with no net magnetic mo- 
ment, or materials with no internal demagnetizing field, 
no correction due to sample boundary is needed, and the 
simple Ewald sum results may be used. This is equiva- 
lent to simulating a long thin "needle" of the bulk ma- 
terial. However, if one wishes to simulate a material in 
which the unit cell has a net magnetic moment or inter- 
nal demagnetizing fields, then we must modify the Ewald 
sum to take into account the necessary boundary effects. 
This is especially important in our simulation because 
the long-range nature of the dipole-dipole interactions 
greatly accentuates these effects. 

The approach described by de Leeuw et al£& is to in- 
clude a boundary term in the Ewald sum of the form 



4tt 



2/2' + 1 J L 3 



(C2) 



B = Mo (H + M) , 



(CI) 



where [ii is the magnetic dipole moment of a spin, and L 
is the system linear dimension. Physically, the inclusion 
of this term corresponds to the consideration of a region 
external to the spherical Ewald boundary (see discussion 
m Ref. Hj). This external region is a continuum with 
magnetic permeability constant ranging from fi' = 1 to 
infinity. Because of the nature of our Ewald sum, the 
fx' = 1 case will in effect simulate the bulk of a spheri- 
cal sample surrounded by a vacuum. The fj/ = oo case 
corresponds to simulating a bulk sample which is needle- 
like and parallel to an applied B, and hence contains no 
internal demagnetizing field. 

In summary, to make a meaningful comparison be- 
tween simulation and experiment within the dipolar spin 
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FIG. 29: Specific heat curves for Dy2Ti2 07 in an applied 
field h//[110], in simulations without (line) and with (circles) 
the boundary term (BT) in the Ewald energy summation. 
These simulations employed single spin flip dynamics in the 
Monte Carlo. The significance of this data in the context 
of experimental measurements on Dy2Ti207 is discussed in 
detail in Section rvil 



ice model in an applied magnetic field, one of two sce- 
narios must happen: 

1. Simulations are performed using the regular Ewald 
summation method, and bulk demagnetization ef- 
fects are accounted for by experimentalists. 

2. Simulations are performed with the inclusion of a 



boundary term Eq. (|C2(1 . Experimentalists are re- 
stricted to measurements on spherical samples to 
make quantitative comparisons. However, mea- 
surements on other sample shapes (with approx- 
imately constant internal fields) may allow some 
qualitative comparison. 

The inclusion of the boundary term is a non-trivial 
matter in many simulations. For example, it will promote 
effects such as domain formation in simulations of global 
Ising ferromagnets. It is therefore always important to 
check ground state configurations of system where the 
term Eq. (|C2(1 is absent against those where is has been 
included, in order not to miss any important secondary 
effects. 

As a means of addressing some of these issues, we 
present some preliminary results on Monte Carlo simula- 
tions of the dipolar spin ice model. Fig. 1291 shows specific 
heat curves for Dy2Ti207 in an applied magnetic field, 
with and without the inclusion of the boundary term, 
Eq. (|C2[) with // = 1. The results presented in Fig. TM 
were performed using single spin flip dynamics on a sys- 
tem of size L = 3. As we see, the boundary term does 
not effect the h = specific heat curve. This is consis- 
tent with the understanding that the boundary term is 
only necessary in Monte Carlo programs where the sim- 
ulation cell has a net magnetic moment (which is not 
true in general for the spin ice manifold). For moderate 
fields, we see that the boundary term significantly alters 
the shape of the specific heat curve, as expected, since 
the simulation cell is acquiring a net magnetic moment. 
For very large fields, the boundary term begins to lose 
its effect, as the field term (Eq. (12911 ) becomes dominant 
in the Hamiltonian. 
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